{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "858c5f0b",
   "metadata": {},
   "source": [
    "# LU Decompostion\n",
    "**강좌**: *수치해석*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57a6f719",
   "metadata": {},
   "source": [
    "## LU 분해 개요\n",
    "- Gauss Elimination은 선형 방정식 $Ax=b$ 를 Upper triangular matrix $U$를 이용한 $Ux=c$ 로 바꿔서 푸는 방법이다. \n",
    "- 이론적으로 $A=LU$ 로 분해해서 푸는 것과 같다.\n",
    "\n",
    "$$\n",
    "    Ax=LUx=b\\\\\n",
    "    Ux=c, Lc=b\n",
    "$$\n",
    "\n",
    "- $A=LU$ 로 분해해놓으면 $b$ 벡터가 달라져도 Substitution 과정만으로 빠르게 계산할 수 있다."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "916ae787",
   "metadata": {},
   "source": [
    "## 이론\n",
    "- Gauss Elimination 과정을 Matrix Operation으로 생각하면 다음과 같다.\n",
    "    * 아래 예제를 생각하자.\n",
    "        \n",
    "        $$\n",
    "        A = \\begin{bmatrix}\n",
    "        2 & 1 & 1 \\\\\n",
    "        4 & 6 & 0 \\\\\n",
    "        -2 & 7 & 2\n",
    "        \\end{bmatrix}\n",
    "        $$\n",
    "        \n",
    "    * $a_{21}$ 을 제거하기 위한 Row Operation ($(2) - 2\\times(1)$) 연산 Matrix\n",
    "    \n",
    "        $$\n",
    "        E_{21}=\\left[\\begin{array}{ccc}\n",
    "        1 & 0 & 0\\\\\n",
    "        -2 & 1 & 0\\\\\n",
    "        0 & 0 & 1\n",
    "        \\end{array}\\right]$$\n",
    "        \n",
    "    * $a_{31}$ 을 제거하기 위한 Row Operation ($(3) + 1\\times(1)$) 연산 Matrix\n",
    "    \n",
    "        $$E_{31}=\\left[\\begin{array}{ccc}\n",
    "        1 & 0 & 0\\\\\n",
    "        0 & 1 & 0\\\\\n",
    "        l & 0 & 1\n",
    "        \\end{array}\\right]$$       \n",
    "        \n",
    "    * $a_{32}$ 을 제거하기 위한 Row Operation ($(2) + 1\\times(1)$) 연산 Matrix\n",
    "    \n",
    "        $$E_{32}=\\left[\\begin{array}{ccc}\n",
    "        1 & 0 & 0\\\\\n",
    "        0 & 1 & 0\\\\\n",
    "        0 & l & 1\n",
    "        \\end{array}\\right]$$       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b51cb2a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "def ematrix(n,i,j,d):\n",
    "    \"\"\"\n",
    "    Elimination matrix\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    n : integer\n",
    "        size of matrix\n",
    "    i : integer\n",
    "        row of pivot\n",
    "    j : intger\n",
    "        row of elimination\n",
    "    d : float\n",
    "        elimination coefficeint\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    mat : array\n",
    "        elimination matrix\n",
    "    \"\"\"\n",
    "    # Make Identity matrix (size n)\n",
    "    mat = np.eye(n)        \n",
    "    \n",
    "    # Assign elimination coefficient\n",
    "    mat[j,i] = d\n",
    "    \n",
    "    return mat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d0e09b1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Problem\n",
    "A = np.array([[2,1,1],[4,-6,0],[-2,7,2]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "495dbf38",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 2.  1.  1.]\n",
      " [ 0. -8. -2.]\n",
      " [-2.  7.  2.]]\n"
     ]
    }
   ],
   "source": [
    "# One step (E21)\n",
    "E21 = ematrix(3, 0, 1, -2)\n",
    "print(E21 @ A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "dc2d9dd7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 2.  1.  1.]\n",
      " [ 0. -8. -2.]\n",
      " [ 0.  0.  1.]]\n"
     ]
    }
   ],
   "source": [
    "# For all steps\n",
    "E31 = ematrix(3, 0, 2, 1)\n",
    "E32 = ematrix(3, 1, 2, 1)\n",
    "U = E32 @ E31 @ E21 @ A \n",
    "print(U)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4ad76630",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.],\n",
       "       [-2.,  0.,  0.]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Not commutative\n",
    "E32 @ E31 @ E21  - E21 @ E31 @ E32"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8ce03a5",
   "metadata": {},
   "source": [
    "- $E_{32} E_{31} E_{21} A = U$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "13483e20",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Undo Gauss Eliminate\n",
      "[[ 2.  1.  1.]\n",
      " [ 4. -6.  0.]\n",
      " [-2.  7.  2.]]\n",
      "Lower Triangular matrix\n",
      "[[ 1.  0.  0.]\n",
      " [ 2.  1.  0.]\n",
      " [-1. -1.  1.]]\n"
     ]
    }
   ],
   "source": [
    "iE32 = ematrix(3, 1, 2, -1)\n",
    "iE31 = ematrix(3, 0, 2, -1)\n",
    "iE21 = ematrix(3, 0, 1, 2)\n",
    "\n",
    "iE = iE21 @ iE31 @ iE32\n",
    "print(\"Undo Gauss Eliminate\")\n",
    "print(iE @ U)\n",
    "\n",
    "print(\"Lower Triangular matrix\")\n",
    "print(iE)\n",
    "\n",
    "L =iE"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c9f9a20",
   "metadata": {},
   "source": [
    "- LU decomposition\n",
    "    * Formula\n",
    "    \n",
    "        $$\n",
    "        LU = \\begin{bmatrix}\n",
    "        u_{11} & u_{12} & u_{13} \\\\\n",
    "        0 & u_{22} & u_{23} \\\\\n",
    "        0 & 0 & u_{33} \\\\\n",
    "        \\end{bmatrix}\n",
    "        \\begin{bmatrix}\n",
    "        1 & 0 & 0 \\\\\n",
    "        l_{21} & 1 & 0 \\\\\n",
    "        l_{31} & l_{23} & 1\\\\\n",
    "        \\end{bmatrix}\n",
    "        $$\n",
    "        \n",
    "        - $l_{ij} = a'_{ij} / a'_{ii}$\n",
    "        \n",
    "    * Save in a matrix\n",
    "        \n",
    "        $$\n",
    "        M = \\begin{bmatrix}\n",
    "        u_{11} & u_{12} & u_{13} \\\\\n",
    "        l_{21} & u_{22} & u_{23} \\\\\n",
    "        l_{31} & l_{23} & u_{33} \\\\\n",
    "        \\end{bmatrix}\n",
    "        $$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e15cfd6",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "- Make $A \\rightarrow M$\n",
    "    \n",
    "    $$\n",
    "    U = \\left[\\begin{array}{ccc}\n",
    "    2 & 1 & 1\\\\\n",
    "    0 & -8 & -2\\\\\n",
    "    0 & 0 & 1\n",
    "    \\end{array}\\right],\n",
    "    L = \\left[\\begin{array}{ccc}\n",
    "    1 & 0 & 0\\\\\n",
    "    2 & 1 & 0\\\\\n",
    "    -1 & -1 & 1\n",
    "    \\end{array}\\right]\n",
    "    $$\n",
    "    \n",
    "    $$\n",
    "    M = \\left[\\begin{array}{ccc}\n",
    "    2 & 1 & 1\\\\\n",
    "    2 & -8 & -2\\\\\n",
    "    -1 & -1 & 1\n",
    "    \\end{array}\\right],\n",
    "    $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "01c8e418",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[2, 1, 1], [4, -6, 0], [-2, 7, 2]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "9f03bdac",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.0 [ 2 -8 -2]\n"
     ]
    }
   ],
   "source": [
    "# First pivot a[0, 0]\n",
    "# eliminate a[1, 0]\n",
    "i = 0\n",
    "j = 1\n",
    "ratio = A[j, i] / A[i, i]\n",
    "\n",
    "A[j, i] = ratio\n",
    "A[j, i+1:] = A[j, i+1:] - ratio*A[i, i+1:]\n",
    "print(ratio, A[j])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "d966d262",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.0 [-1  8  3]\n"
     ]
    }
   ],
   "source": [
    "# eliminate a[2, 0]\n",
    "j = 2\n",
    "ratio = A[j, i] / A[i, i]\n",
    "\n",
    "A[j, i] = ratio\n",
    "A[j, i+1:] = A[j, i+1:] - ratio*A[i, i+1:]\n",
    "print(ratio, A[j])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "64d8ce3c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.0 [-1 -1  1]\n"
     ]
    }
   ],
   "source": [
    "# next pivot a_{2,2}\n",
    "# eliminate a_{3, 2}\n",
    "i = 1\n",
    "j = 2\n",
    "\n",
    "ratio = A[j, i] / A[i, i]\n",
    "\n",
    "A[j, i] = ratio\n",
    "A[j, i+1:] = A[j, i+1:] - ratio*A[i, i+1:]\n",
    "print(ratio, A[j])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "bb1ba558",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 2  1  1]\n",
      " [ 2 -8 -2]\n",
      " [-1 -1  1]]\n"
     ]
    }
   ],
   "source": [
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e1a3511",
   "metadata": {},
   "source": [
    "- Substibution\n",
    "    * Forward : $Lc=b$\n",
    "    * Backward: $Ux=c$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "f64d293f",
   "metadata": {},
   "outputs": [],
   "source": [
    "b = np.array([5, -2, 9])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "13ddb204",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Forward (b to c)\n",
    "# First row\n",
    "b[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "8fa59f5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Second row\n",
    "b[1] = b[1] - A[1, 0]*b[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "319babde",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  5 -12   2]\n"
     ]
    }
   ],
   "source": [
    "# Third row\n",
    "i = 2\n",
    "\n",
    "for j in range(i):\n",
    "    b[i] -= A[i, j]*b[j]\n",
    "    \n",
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "f53ab74d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Back substitution (Same as Gauss Elimination)\n",
    "x = np.empty_like(b)\n",
    "\n",
    "# Third row\n",
    "i = 2\n",
    "x[i] = b[i] / A[i,i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "16509053",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Second row\n",
    "i = 1\n",
    "x[i] = (b[i] - A[i, i+1]*x[i+1]) / A[i, i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "1a8a1a06",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 1 2]\n"
     ]
    }
   ],
   "source": [
    "# First row\n",
    "n = 3\n",
    "i = 0\n",
    "x[i] = b[i]\n",
    "\n",
    "for j in range(i+1, n):\n",
    "    x[i] -= A[i, j]*x[j]\n",
    "    \n",
    "x[i] /= A[i,i]\n",
    "\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "024c240c",
   "metadata": {},
   "source": [
    "### DIY\n",
    "다음 2개의 함수를 만드시오\n",
    "* LU 분해: LUdecomp\n",
    "* Substitution : LUsubs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9700a890",
   "metadata": {},
   "source": [
    "## 역행렬\n",
    "역행렬은 $Ax_i = e_i$ 를 푸는 과정이다.\n",
    "\n",
    "$$\n",
    "A X = I = A \n",
    "\\begin{bmatrix}\n",
    "x_1 & x_2 & ... & x_n\n",
    "\\end{bmatrix} =\n",
    "\\begin{bmatrix}\n",
    "e_1 & e_2 & ... & e_n\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "### 예제\n",
    "다음 $A$ 행렬의 역행렬을 구하시오.\n",
    "\n",
    "$$\n",
    "A = \\begin{bmatrix}\n",
    "2 & 1 & 1 \\\\\n",
    "4 & 6 & 0 \\\\\n",
    "-2 & 7 & 2\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "        \n",
    "LU 분해 결과는 아래와 같다.\n",
    "\n",
    "$$\n",
    "M = \\left[\\begin{array}{ccc}\n",
    "2 & 1 & 1\\\\\n",
    "2 & -8 & -2\\\\\n",
    "-1 & -1 & 1\n",
    "\\end{array}\\right],\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "1ed5aa61",
   "metadata": {},
   "outputs": [],
   "source": [
    "# LU 분해 결과 (이전 계산)\n",
    "A = np.array([[ 2,  1,  1],  [ 2, -8, -2],  [-1, -1,  1]], dtype=np.float64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "c7c076ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Allocate Identity and inverse Matrix\n",
    "I = np.eye(3)\n",
    "Ainv = np.empty_like(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "72e9f0df",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.75  0.5  -1.  ]\n"
     ]
    }
   ],
   "source": [
    "# For first column\n",
    "b = I[:, 0].copy()\n",
    "x = Ainv[:, 0]\n",
    "\n",
    "# Forward (b to c)\n",
    "# Second row\n",
    "b[1] = b[1] - A[1, 0]*b[0]\n",
    "\n",
    "# Third row\n",
    "i = 2\n",
    "\n",
    "for j in range(i):\n",
    "    b[i] -= A[i, j]*b[j]\n",
    "    \n",
    "# Back substitution (Same as Gauss Elimination)\n",
    "# Third row\n",
    "i = 2\n",
    "x[i] = b[i] / A[i,i]\n",
    "\n",
    "# Second row\n",
    "i = 1\n",
    "x[i] = (b[i] - A[i, i+1]*x[i+1]) / A[i, i]\n",
    "\n",
    "# First row\n",
    "n = 3\n",
    "i = 0\n",
    "x[i] = b[i]\n",
    "\n",
    "for j in range(i+1, n):\n",
    "    x[i] -= A[i, j]*x[j]\n",
    "    \n",
    "x[i] /= A[i,i]\n",
    "\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "1eb645c5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.3125 -0.375   1.    ]\n"
     ]
    }
   ],
   "source": [
    "# For first column\n",
    "b = I[:, 1].copy()\n",
    "x = Ainv[:, 1]\n",
    "\n",
    "# Forward (b to c)\n",
    "# Second row\n",
    "b[1] = b[1] - A[1, 0]*b[0]\n",
    "\n",
    "# Third row\n",
    "i = 2\n",
    "\n",
    "for j in range(i):\n",
    "    b[i] -= A[i, j]*b[j]\n",
    "    \n",
    "# Back substitution (Same as Gauss Elimination)\n",
    "# Third row\n",
    "i = 2\n",
    "x[i] = b[i] / A[i,i]\n",
    "\n",
    "# Second row\n",
    "i = 1\n",
    "x[i] = (b[i] - A[i, i+1]*x[i+1]) / A[i, i]\n",
    "\n",
    "# First row\n",
    "n = 3\n",
    "i = 0\n",
    "x[i] = b[i]\n",
    "\n",
    "for j in range(i+1, n):\n",
    "    x[i] -= A[i, j]*x[j]\n",
    "    \n",
    "x[i] /= A[i,i]\n",
    "\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "940f54b2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.375 -0.25   1.   ]\n"
     ]
    }
   ],
   "source": [
    "# For first column\n",
    "b = I[:, 2].copy()\n",
    "x = Ainv[:, 2]\n",
    "\n",
    "# Forward (b to c)\n",
    "# Second row\n",
    "b[1] = b[1] - A[1, 0]*b[0]\n",
    "\n",
    "# Third row\n",
    "i = 2\n",
    "\n",
    "for j in range(i):\n",
    "    b[i] -= A[i, j]*b[j]\n",
    "    \n",
    "# Back substitution (Same as Gauss Elimination)\n",
    "# Third row\n",
    "i = 2\n",
    "x[i] = b[i] / A[i,i]\n",
    "\n",
    "# Second row\n",
    "i = 1\n",
    "x[i] = (b[i] - A[i, i+1]*x[i+1]) / A[i, i]\n",
    "\n",
    "# First row\n",
    "n = 3\n",
    "i = 0\n",
    "x[i] = b[i]\n",
    "\n",
    "for j in range(i+1, n):\n",
    "    x[i] -= A[i, j]*x[j]\n",
    "    \n",
    "x[i] /= A[i,i]\n",
    "\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "8ea60888",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.75   -0.3125 -0.375 ]\n",
      " [ 0.5    -0.375  -0.25  ]\n",
      " [-1.      1.      1.    ]]\n"
     ]
    }
   ],
   "source": [
    "print(Ainv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "75c70c63",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0., 0.],\n",
       "       [0., 1., 0.],\n",
       "       [0., 0., 1.]])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Validate\n",
    "A = np.array([[2,1,1],[4,-6,0],[-2,7,2]])\n",
    "A @ Ainv"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "895b7896",
   "metadata": {},
   "source": [
    "### DIY\n",
    "앞서 만든 LU 분해 함수를 이용하여 역행렬 계산 함수를 만드시오"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84e067bf",
   "metadata": {},
   "source": [
    "## 특수 행렬 분해\n",
    "### Cholesky Decomposition\n",
    "- $LU$ 분해에서 $U$ 행렬은 Diagonal 항과 Off-diagonal 항으로 분리할 수 있다.\n",
    "\n",
    "    $$\n",
    "        A = LU = LDU\n",
    "    $$\n",
    "    \n",
    "- $A$ 행렬이 대칭 (symmetric) 일 경우 다음 관계를 만족한다.\n",
    "\n",
    "    $$\n",
    "        A = A^T \\\\\n",
    "        LDU = (LDU)^T = U^T D L^T\n",
    "    $$\n",
    "    \n",
    "    * 즉 $L = U^T$ 를 만족한다.\n",
    "    * $A=LDL^T$\n",
    "    \n",
    "- Cholesky 분해법\n",
    "    * A 가 Positive Definite 일 때 (symmetric, all positive pivots)\n",
    "    \n",
    "    $$\n",
    "        A = LDL^T = (LD^{1/2}) (D^{1/2} L^T) = (LD^{1/2}) (LD^{1/2})^T\n",
    "    $$\n",
    "    \n",
    "- 구현\n",
    "    \n",
    "    $$\n",
    "    \\begin{align}\n",
    "    l_{ki} &= \\frac{a_{ki} - \\sum_{j=1}^{i-1} l_{ij} l_{kj}}{l_{ii}} \\textrm{ for } i=1,2,...k-1 \\\\\n",
    "    l_{kk} &= \\sqrt{a_{kk} - \\sum_{j=1}^{k-1} l_{kj}^2}\n",
    "    \\end{align}\n",
    "    $$\n",
    "    \n",
    "#### 예제\n",
    "다음 대칭 행렬을 Cholesky 분해하시오.\n",
    "\n",
    "$$\n",
    "A = \\begin{bmatrix}\n",
    "6 & 15 & 55 \\\\\n",
    "15 & 55 & 225 \\\\\n",
    "55 & 225 & 979\n",
    "\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "dab369ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[6, 15, 55], [15, 55, 225], [55, 225, 979]], dtype=np.float64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "244f0142",
   "metadata": {},
   "outputs": [],
   "source": [
    "L = np.zeros_like(A)\n",
    "\n",
    "# first row\n",
    "k = 0\n",
    "L[k, k] = np.sqrt(A[k, k])\n",
    "\n",
    "# Second row\n",
    "k = 1\n",
    "i = 0\n",
    "L[k, i] = A[k, i] / L[i, i]\n",
    "L[k, k] = np.sqrt(A[k, k] - L[k, i]**2)\n",
    "\n",
    "# Third row\n",
    "k = 2\n",
    "i = 0 \n",
    "L[k, i] = A[k, i] / L[i, i]\n",
    "\n",
    "i = 1\n",
    "L[k, i] = (A[k, i] - sum(L[i, j]*L[k, j] for j in range(i)))/L[i,i]\n",
    "\n",
    "L[k ,k] = np.sqrt(A[k, k] - sum(L[k, j]**2 for j in range(k)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "7eebca22",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 2.44948974,  0.        ,  0.        ],\n",
       "       [ 6.12372436,  4.18330013,  0.        ],\n",
       "       [22.45365598, 20.91650066,  6.11010093]])"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "032d669a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  6.,  15.,  55.],\n",
       "       [ 15.,  55., 225.],\n",
       "       [ 55., 225., 979.]])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Validate\n",
    "L @ L.T"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6efc5570",
   "metadata": {},
   "source": [
    "- LU decomposition과 비교"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "72cf5f99",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 6.         15.         55.        ]\n",
      " [ 2.5        17.5        87.5       ]\n",
      " [ 9.16666667  5.         37.33333333]]\n"
     ]
    }
   ],
   "source": [
    "# LU decomposition\n",
    "M = A.copy()\n",
    "i = 0\n",
    "j = 1\n",
    "ratio = M[j, i] / M[i, i]\n",
    "\n",
    "M[j, i] = ratio\n",
    "M[j, i+1:] = A[j, i+1:] - ratio*M[i, i+1:]\n",
    "\n",
    "j = 2\n",
    "ratio = M[j, i] / M[i, i]\n",
    "\n",
    "M[j, i] = ratio\n",
    "M[j, i+1:] = M[j, i+1:] - ratio*M[i, i+1:]\n",
    "\n",
    "i = 1\n",
    "j = 2\n",
    "ratio = M[j, i] / M[i, i]\n",
    "\n",
    "M[j, i] = rati\n",
    "M[j, i+1:] = M[j, i+1:] - ratio*M[i, i+1:]\n",
    "\n",
    "print(M)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "112200aa",
   "metadata": {},
   "source": [
    "- $A = LU=LDU'$\n",
    "$$\n",
    "L = \\begin{bmatrix}\n",
    "1 & 0 & 0 \\\\\n",
    "2.5 & 1 & 0 \\\\\n",
    "9.1667 & 5 & 1\n",
    "\\end{bmatrix}\n",
    "U = \\begin{bmatrix}\n",
    "6 & 1.5 & 55 \\\\\n",
    "0 & 17.5 & 87.5 \\\\\n",
    "0 & 0 & 37.3333\n",
    "\\end{bmatrix}\n",
    "= \\begin{bmatrix}\n",
    "6 & 0 & 0 \\\\\n",
    "0 & 17.5 & 0 \\\\\n",
    "0 & 0& 37.3333\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "1 & 2.5 & 9.1667 \\\\\n",
    "0 & 1 &  5 \\\\\n",
    "0 & 0 & 1\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "- $L\\sqrt{D}$\n",
    "$$\n",
    "L \\sqrt{D} = \\begin{bmatrix}\n",
    "\\sqrt{6} & 0 & 0 \\\\\n",
    "2.5 \\sqrt{6} & \\sqrt{17.5} & 0 \\\\\n",
    "9.1667 \\sqrt{6} & 5 \\sqrt{17.5} & \\sqrt{37.3333}\n",
    "\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c509bc0",
   "metadata": {},
   "source": [
    "### 삼중대각 행렬\n",
    "- 삼중 대각 행렬 (Tri-diagonal matrix)는 Diagonal과 그 바로 옆에 1개씩만 값이 있는, 즉 3의 대역폭을 갖는 행렬이다.\n",
    "\n",
    "- 실제 편미분 방정식 해석에서 자주 보이는 행렬이다.\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "b_1 & c_1 & & & & & \\\\\n",
    "a_2 & b_2 & c_2 & & & & \\\\\n",
    "& a_3 & b_3 & c_3 & & & \\\\\n",
    "& & & & & & & \\\\\n",
    "& & & & a_{n-1} & b_{n-1} & c_{n-1} \\\\\n",
    "& & & & & a_{n} & b_{n}\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "x_1 \\\\ x_2 \\\\ x_3 \\\\ \\\\ x_{n-1} \\\\ x_n\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "d_1 \\\\ d_2 \\\\ d_3 \\\\ \\\\ d_{n-1} \\\\ d_n\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "- (LU Decomposition + Forward substitution), Backward substitution 으로 효과적으로 계산할 수 있다.\n",
    "\n",
    "- Algorithm\n",
    "    * For $i = 2, 3,...,n$\n",
    "        \n",
    "        $$\n",
    "        \\begin{align}\n",
    "            w &= \\frac{a_i}{b_{i-1}} \\\\\n",
    "            b_i &= b_i - w c_{i-1} \\\\\n",
    "            d_i &= d_i - w d_{i-1}\n",
    "         \\end{align}\n",
    "        $$\n",
    "        \n",
    "    * Back substitution\n",
    "    \n",
    "        $$\n",
    "        \\begin{align}\n",
    "        x_n &= \\frac{d_n}{b_n} \\\\\n",
    "        x_i & = \\frac{d_i - c_i x_{i+1}}{b_i} \\textrm{ for } i=n-1, n-2,...,1        \n",
    "        \\end{align}\n",
    "        $$\n",
    "\n",
    "#### 예제\n",
    "아래 삼중대각 행렬의 해를 구하시오.\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "2.04 & -1 & & \\\\\n",
    "-1 & 2.04 & -1 & \\\\\n",
    "& -1 & 2.04 & -1 \\\\\n",
    "& & -1 & 2.04\n",
    "\\end{bmatrix}\n",
    "T \n",
    "= \n",
    "\\begin{bmatrix}\n",
    "40.8 \\\\ 0.8 \\\\ 0.8 \\\\ 200.8\n",
    "\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "f7df4b3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[2.04, -1, 0, 0], [-1, 2.04, -1, 0], [0, -1, 2.04, -1], [0, 0, -1, 2.04]])\n",
    "d = np.array([40.8, 0.8, 0.8, 200.8])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "fa9a4e0c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 2.04, -1.  ,  0.  ,  0.  ],\n",
       "       [-1.  ,  2.04, -1.  ,  0.  ],\n",
       "       [ 0.  , -1.  ,  2.04, -1.  ],\n",
       "       [ 0.  ,  0.  , -1.  ,  2.04]])"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "1f0ebcf3",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([0., -1., -1. , -1.])\n",
    "b = np.array([2.04, 2.04, 2.04, 2.04])\n",
    "c = np.array([-1., -1., -1., 0.])\n",
    "\n",
    "n = len(a)\n",
    "x = np.empty_like(a)\n",
    "\n",
    "# Forward sweep\n",
    "for i in range(1, n):\n",
    "    w = a[i] / b[i-1]\n",
    "    b[i] -= w*c[i-1]\n",
    "    d[i] -= w*d[i-1]\n",
    "\n",
    "# Backward sweep    \n",
    "x[-1] = d[-1] / b[-1]\n",
    "for i in range(n-2, -1, -1):\n",
    "    x[i] = (d[i] - c[i]*x[i+1]) / b[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "a5150aa6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 65.96983437,  93.77846211, 124.53822833, 159.47952369])"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "5db24d68",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 40.8,   0.8,   0.8, 200.8])"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Validation\n",
    "A @ x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a39720d1",
   "metadata": {},
   "source": [
    "### DIY\n",
    "- Cholesky 분해 계산 함수를 만드시오.\n",
    "- 삼중 대각 행렬 계산 함수를 만드시오. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd25f5b8",
   "metadata": {},
   "source": [
    "## Scipy 활용\n",
    "\n",
    "`scipy.linalg` 모듈은 다양한 행렬 계산 알고리즘을 제공함\n",
    "* `numpy.linalg` 는 제한된 기능을 제공\n",
    "\n",
    "**참고**\n",
    "* https://scipy-lectures.org/intro/scipy.html#linear-algebra-operations-scipy-linalg\n",
    "* https://docs.scipy.org/doc/scipy/tutorial/linalg.html\n",
    "\n",
    "### Solve linear system\n",
    "\n",
    "- `linalg.solve` : $Ax=b$ 해석\n",
    "- `linalg.inv` : $A^{-1}$ 계산\n",
    "- `linalg.det` : $det(A)$ 계산\n",
    "- `linalg.norm` : 벡터, 행렬의 Norm 계산\n",
    "\n",
    "### Decomposition\n",
    "- LU (Lower Upper) Decomposition\n",
    "- Singual Value Decomposition (SVD)\n",
    "- QR Decomposition\n",
    "\n",
    "### Eigenvalue \n",
    "- `linalg.eig` : Eigenvalue, Eigenvector 계산\n",
    "    * $Ax = \\lambda x$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "fbecfb26",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.linalg.solve?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c1f5088",
   "metadata": {},
   "source": [
    "### 예제\n",
    "양 끝단이 고정된 보의 방정식은 다음과 같다.\n",
    "\n",
    "$$\n",
    "u''(x) = 1, u(0) = 0, u(1) = 0.\n",
    "$$\n",
    "\n",
    "이 미분방정식의 이론해는 다음과 같다.\n",
    "\n",
    "$$\n",
    "u(x) = -\\frac{1}{2}x + \\frac{1}{2}x^2.\n",
    "$$\n",
    "\n",
    "수치적으로 이를 해석하는 경우 $[0,1]$ 구간을 *n* 등분하고, 각 지점의 값을 $u_i = u(x_i)$ 라 하면 위 미분 방정식으 다음 형태로 바뀐다.\n",
    "\n",
    "$$\\frac{1}{\\Delta x^2}\n",
    "\\begin{bmatrix}\n",
    "-2 & 1 & 0 &... &0 \\\\\n",
    "1 & -2& 1 & 0 ... & 0 \\\\\n",
    "0 & 1 & -2 & 1 ... & 0 \\\\\n",
    "... & ... & ... & ... & ... \\\\\n",
    "0 & 0 & 1 & -2 & 1 \\\\\n",
    "0 & 0 & 0 & 1 & -2\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "u_1 \\\\\n",
    "u_2 \\\\\n",
    "u_3 \\\\\n",
    "... \\\\\n",
    "u_{n-2} \\\\\n",
    "u_{n-1}\n",
    "\\end{bmatrix}\n",
    "= \\begin{bmatrix}\n",
    "f(x_1) \\\\\n",
    "f(x_2) \\\\\n",
    "f(x_3) \\\\\n",
    "... \\\\\n",
    "f(x_{n-2})\\\\\n",
    "f(x_{n-1})\\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "이 문제를 수치적으로 해석해보자."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "c29bcc98",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of division\n",
    "n = 51\n",
    "x = np.linspace(0, 1, n+1)\n",
    "h = np.diff(x)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "ba41eba7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Matrix system\n",
    "A = np.diag([-2]*(n-1)) + np.diag([1]*(n-2), -1) + np.diag([1]*(n-2), 1)\n",
    "b = np.array([1]*(n-1))*h**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "8e6592d7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[-2,  0,  0,  0],\n",
       "        [ 0, -2,  0,  0],\n",
       "        [ 0,  0, -2,  0],\n",
       "        [ 0,  0,  0, -2]]),\n",
       " array([[0, 0, 0, 0],\n",
       "        [1, 0, 0, 0],\n",
       "        [0, 1, 0, 0],\n",
       "        [0, 0, 1, 0]]),\n",
       " array([[0, 1, 0, 0],\n",
       "        [0, 0, 1, 0],\n",
       "        [0, 0, 0, 1],\n",
       "        [0, 0, 0, 0]]))"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.diag([-2]*(n-1)), np.diag([1]*(n-2), -1), np.diag([1]*(n-2), 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "ed8279c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solve\n",
    "#u = np.linalg.solve(A, b)\n",
    "from scipy import linalg\n",
    "u = linalg.solve(A, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "41d10132",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f4267518410>"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1QAAAJqCAYAAAAsZgE0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAABcSAAAXEgFnn9JSAACjTklEQVR4nOzdd3Qc1d0+8Gd7Ue9dsmVZbpKLLFtyl3HFNqaGEOIYkpDkJS8BB5KQkIB/kAKEYHhJ4lQSsKkGbNx7k7tsuclNrrJs9d5XW39/LJqVULH63V09n3M4Z+7szOhZeyz2u/fOvTKbzWYDERERERERdZlcdAAiIiIiIiJXxYKKiIiIiIiom1hQERERERERdRMLKiIiIiIiom5iQUVERERERNRNLKiIiIiIiIi6iQUVERERERFRN7GgIiIiIiIi6iYWVERERERERN3EgoqIiIiIiKibWFARERERERF1EwsqIiIiIiKibmJBRURERERE1E0sqIiIiIiIiLpJKTrAQFVYWAibzSbs5wcGBgIASktLhWUg58X7g9rDe4Paw3uD2sN7gzoi+v6QyWQIDQ3t0TVYUAlis9mEFlTNcxC1h/cHtYf3BrWH9wa1h/cGdcSV7w8O+SMiIiIiIuomFlRERERERETdxIKKiIiIiIiom1hQERERERERdRMLKiIiIiIiom5iQUVERERERNRNLKiIiIiIiIi6iQUVERERERFRN7GgIiIiIiIi6iYWVERERERERN3EgoqIiIiIiKibWFARERERERF1EwsqIiIiIiKibmJBRURERERE1E0sqIiIiIiIiLpJKTrA1xmNRnz55Zc4dOgQSktL4enpiTFjxuCb3/wmAgICunSturo6fPbZZ8jIyEBlZSV8fX0xYcIEPPzww/Dw8GjzHKvViq1bt2LPnj0oLCyEVqvFqFGj8PDDDyMyMrI33iIREREREbkJp+qhMhqN+O1vf4vPP/8cBoMBycnJCAgIwL59+/D888+jsLCw09eqqanBCy+8gC1btkChUGDChAnQ6XTYunUrfvWrX6GmpqbVOTabDW+//Tbef/99lJeXIykpCVFRUTh27Bh++ctf4sqVK735domIiIiIyMU5VQ/VunXrkJ2djfj4ePzmN7+BVqsFAGzatAmrVq3C3/72N7z88sudutb777+PgoICTJw4ET/96U+hUCgAAP/5z3+wbds2vP/++3jqqadanLN3714cPXoUYWFhePnll+Hr6wsAOHr0KFasWIF33nkHb7/9tnQtIiIiIiIa2Jymh8psNmPbtm0AgO9///tSMQUAixYtQkxMDC5evIjr16/f8VqVlZU4cOAAFAoFnnjiiRYF0He+8x14e3vj4MGDqKysbHHepk2bAADf/va3pWIKAFJTU5GcnIyioiIcP368B++SiIiIiIjcidMUVJcuXUJdXR1CQkIwePDgVq+npKQAAE6cOHHHa506dQo2mw0jR45sURgBgEqlwvjx42G1WnH69Glpf3FxMW7fvg21Wo2kpKRW10xNTQUAZGZmduFdERERERGRO3OagurmzZsA0GYxBQCxsbEtjuvJtZr25+TkSPuatqOioqBUth4J2XROZ36+s7PZbLhWWod/Hs5BZl6t6DhERERENACl51Tjv8dycauiQXSUHnGaZ6hKS0sBoN2Z/Pz9/Vsc15lrNZ3zdU0/o/m17vTz2zqnI88++2yrfWq1Gq+99hoAIDAwsFPX6Qt/P5SDVcezAQDThwRg/ti2C08auJq+VAgKChKchJwN7w1qD+8Nag/vDWrP9j15OFdQg38duYmnp8fikaQI0ZG6xWl6qAwGAwBAo9G0+XrTM1VNx/XkWk37GxsbW52jVqs7PKczP9/ZjY3wkbaP5pSjrtEsMA0RERERDTSF1QacK3DMuj0u0qeDo52b0/RQ2Wy2Hr3e1rEymazLObpzTltWrFjR4eulpaVdek+9KUZng49WiSqDGUaLDVvO5CBtsOvexNT7mr5FLCkpEZyEnA3vDWoP7w1qD+8NasvGi2XSdqSvFn6oR0lJ/w/9k8lkCAsL69E1nKaHSqfTAWjZa9Rc0/7ms//d6Vrt9SY1Xat5D1bTdXvj5zs7pVyGGXGOIYcHb1YLTENEREREA83Bm47eqVnxQb3WqSGC0xRUTc8UlZWVtfl6eXl5i+M6c62mc76u6Wc0v9adfn5b57iy2fGO93GqoA61jRaBaYiIiIhooCisMeJKmaPjY3a8az9f5zQFVUxMDADgxo0bbb7etP5U03E9uVbT/ubXGjRoEADg1q1bMJtbP1PUdE50dPQdf74rGBvpCz+9CgBgtgJHb9fc4QwiIiIiop47lOv43DnIX4fYAL3AND3nNAXV8OHDodfrUVRU1GYhdOzYMQBoc42orxs7dixkMhkuXryIqqqqFq+ZTCZkZmZCJpNh3Lhx0v7g4GBERETAaDTi5MmTra559OhRAMD48eO79L6clVIuQ1qzYX+HbrKgIiIiIqK+1/xxk7uGuvZwP8CJCiqlUon58+cDAP7zn/+0eP5p06ZNuHnzJoYPH464uDhp/7Zt27Bs2TJ89NFHLa7l5+eHKVOmwGw249///jcsFsdwtg8++ADV1dWYOnVqq0V/Fy1aBAD48MMPWxRix44dw4kTJxAcHIwJEyb02nsWrfmwvzOFdajmsD8iIiIi6kP51UZcr3DMWTB7mGsP9wOcaJY/AHjggQeQlZWF7OxsPPPMMxg+fDhKS0tx5coVeHl54cc//nGL46urq5Gfn4+KiopW13r88cdx5coVHDt2DMuWLcOQIUNw69Yt3Lp1CyEhIXjsscdanTNz5kycOnUKGRkZWLZsGRITE1FTU4MLFy5ApVLhJz/5SZuL/rqq0eE+8NMpUdFghsUGHL1Vg7lxvqJjEREREZGbat47FeOrwSB/1x7uBzhRDxVgXwNq+fLlePDBB6FWq3H8+HEUFxdjxowZeP311xEaGtrpa3l7e+PVV1/F/PnzYTabkZGRgfr6esyfPx+vvvoqvL29W50jl8vx7LPPYunSpfD390dmZiZyc3MxYcIEvPbaaxg2bFhvvl3hFHIZpkR7Se0DnO2PiIiIiPrQwWbPT02N8ergSNchs4laDGmAKygoELYOFeBYEyL9/E38cmcuAEAuA/57fxx8de7TC0fdwzVDqD28N6g9vDeoPbw3qMmtqkY8tckxV8Lf7onFmCERAMTdH261DhWJMSxIhwC9vYCy2oAjtzg5BRERERH1vubD/WL9NAj3VgtM03tYUA1wcpkMU5sN++Miv0RERETU22w2W4vFfKfGtH78xlWxoKIWN/T54gaU1ZsEpiEiIiIid3OzshG3q41S212enwJYUBGAoQFaBHvYF/m1ATicy2F/RERERNR7mvdODQ3QIsTTPYb7ASyoCPaH8Zp/S3CIBRURERER9RKbzYaDuY7HStypdwpgQUVfaT7s72JJA0rqOOyPiIiIiHruRkUjCmocny2nRLvP81MACyr6SqyfBmFeKqnNYX9ERERE1Buar3U6PFCHIA9VB0e7HhZUBMA+7K/5twVc5JeIiIiIeqr17H7uNdwPYEFFzUxrdoNfKTOgqNbYwdFERERERB27UmZA8VePksgATI5mQUVuLMZXg8hmC6w1/zaBiIiIiKirmk92NjJYhwC9ew33A1hQUTNfn+2Pi/wSERERUXdZbbYWj5G422QUTVhQUQvNZ/u7XtGI29WNAtMQERERkau6WNKAsnozAEAuA6a44fNTAAsq+pooHw0G+2mk9sEcDvsjIiIioq47kOPonRodooevVikwTd9hQUWtTGvWS5V+sxo2m01gGiIiIiJyNRarrcXzU9MGuedwP4AFFbWh+XNUedVG3KjgsD8iIiIi6rwzhXWobrQAAJRyGVKj3HO4H8CCitoQ4qnGsECd1OaaVERERETUFQeazRY9PtwDnmqFwDR9iwUVtan5mlQHcjjsj4iIiIg6x2ix4uitZsP9Ytx3uB/AgoraMTXGG3KZfbuk3oxLpQ1iAxERERGRSziZX4d6kxUAoFHIMCHSU3CivsWCitrkp1MiIVgvtQ9wkV8iIiIi6oT0ZrP7pUR6Qat075LDvd8d9Ujz2VgO3ayGxcphf0RERETUvgaTFcfzaqX21EHuOxlFExZU1K5JUV5o+kKh0mDBueJ6sYGIiIiIyKll3K6B0WL/Et5DLUdSmIfgRH2PBRW1y0ujwLhm/wiad98SEREREX1d89mhJ0V5QaVw/3LD/d8h9UjzWVmO3KqBycJhf0RERETUWk2jBacK6qS2u8/u14QFFXVoYqQX1Ar7dH91RitOFdTe4QwiIiIiGoiO3KqB2T65H3y1CiSG6Ds+wU2woKIO6VRyTIhwTHXJ2f6IiIiIqC3Nh/tNifaComkNHjfHgoruqPlsfxm3a9DY9NUDERERERGAigYzzhU5JjBr/vnR3bGgojsaH+4Bvcp+qxjMNmTcZi8VERERETkcvFmNphV2gvRKDAvUiQ3Uj1hQ0R2pFXKkRDRb5Df9lMA0RERERORsDhw9L21PjfaEXDYwhvsBLKiok6bLSqTtTFkgaqs4hToRERERAYW3CpCtCpLa0zQD63MiCyrqlMTkBHib7NNgmuVKHDmcJTgRERERETmDg8cvSdsRjeUYPDJOYJr+x4KKOkWlVmGyqlJqpxcYxYUhIiIiIqdgs9mQXqmS2tO8GyGXD6wSY2C9W+qRGaMipO1z6hCUFRQLTENEREREouVcuoabmkCpPX38wOqdAlhQURcMGx2PIGMVAMAqk+PgsQuCExERERGRSOmnb0rbcY0liBgcJTCNGCyoqNMUcjmmezjWF0gvGziztxARERFRSxaLGQcMnlJ7elAHB7sxFlTUJdOTYqXtq9oQ5F3JEReGiIiIiITJPnkBJWofAIDcZsXUlJGCE4nBgoq6ZFBcDGKMZVI7/eRVgWmIiIiISJT9F4uk7QRTMQKCAwSmEYcFFXXZdH+rtL2/Tg+rxSIwDRERERH1N1NjIw5ZHQXUjEidwDRisaCiLps2cYS0XaDxx7Wzlzo4moiIiIjczemjZ1Gj0gMAVFYTUiclCk4kDgsq6rKQsECMNDqmTN9/Pk9gGiIiIiLqb/uvV0jbybZSeHrqBaYRiwUVdcu0MMcCbgdNfjAbudAvERER0UDQUFWNDEWo1J4R6ycwjXgsqKhbpqQmQGG1PztVofbCuYyzghMRERERUX84dvgMGhVqAIDebEDSxFGCE4nFgoq6xcfXC+OsJVJ7/9WyDo4mIiIiIneRnt8obU9WVkCjVnVwtPtjQUXdNmOQt7R9VBaCxto6gWmIiIiIqK9VFRThtDpMak8fGdbB0QMDCyrqtgkpo6C12J+dqldqceLwabGBiIiIiKhPHTp6Hha5AgDgZ6rFqNHxghOJx4KKuk2n1SBFUS6102+xh4qIiIjInaU3e8pjmr4OSgXLCf4JUI9MHxYsbWeqw1BTUtLB0URERETkqoquXsdFXbjUnj4uVmAa58GCinpkzLjh8DbXAwBMchWOHD4nOBERERER9YX0E1el7QhTJeLiIgWmcR4sqKhHVAo5pmhrpPaBYovANERERETUF2xWKw7UaqX2NF8zZDKZwETOgwUV9diM0dHSdpYuHKU5NwWmISIiIqLelnPmAm7qHI96zJg4TGAa58KCinps+PBBCDFVAwBsMjkOZFwWnIiIiIiIetP+c7el7aGmUoSHBwlM41xYUFGPyWQyTPcxSu191RrYbDaBiYiIiIiot5gbjdhv9pfaM0IH9kK+X8eCinpFWvJQaTtHF4ycrEsC0xARERFRbzmfcQblam8AgNxmwbRJIwUnci4sqKhXREaFIM5YKrX3Z+UKTENEREREvWXvVce6o+MspfD18RKYxvmwoKJeMyNUKW3vN/rBYjIJTENEREREPWWorsFReYjUThvsLTCNc2JBRb1mWupIyG32adPL1d44n5ElOBERERER9UTGkTNoUNqnS9dZGjFx4ijBiZwPCyrqNX5+3hhrKZHae6+UdnA0ERERETm7fbcN0vYkRTm0WrXANM6JBRX1qrQYx5jaI7JgGOpqBaYhIiIiou6qLCzBKU2Y1J4xPFRgGufFgop6VUpKArSWRgBAg1KL44fOCk5ERERERN1x8Oh5WGUKAIC/qQYJY7mYb1tYUFGv0uo0mCQvk9r7bjcITENERERE3bWvTCZtT9fXQalg6dAW/qlQr0sb5pgJ5pQqFJXFZR0cTURERETO5vaVHFzRNpvdb9xggWmcGwsq6nUJ44bD31QDALDIFTh49JzgRERERETUFftOXpO2YxrLMHhojMA0zo0FFfU6pVKBaboaqb2vpIODiYiIiMipWC0WpNd7SO0Z/haBaZwfCyrqE2ljHd3CV7QhyLueKzANEREREXVW9plsFKl9AQAymxXTU0aIDeTkWFBRnxgcPwjRjc0mpzhxVWAaIiIiIuqsvRfype0EYzGCwoIEpnF+StEBmsvOzsbatWtx+fJlmM1mREZGYt68eUhLS+vW9TIzM7Fhwwbk5OQAAAYNGoTFixdj/PjxrY4tKSnBiRMncPr0aeTl5aG8vBw6nQ6xsbGYN28ekpOTe/DOBh6ZTIY0PxNW1dvb++v0+JbVCrmcNTwRERGRszI2GnHIEiBVCWnhGrGBXIDTFFQZGRlYsWIFbDYbRowYAS8vL5w7dw4rV67EzZs38dhjj3Xpelu2bMF7770HhUKBxMREKJVKnD17Fq+//joef/xxLFiwoMXx77zzDrKzs6FWqzF06FDExcWhqKgIZ86cwZkzZ7Bw4cIuZxjopk0cjtV7y2GTyVGk9kX22csYMXa46FhERERE1I6Tx7JQq7Q/P6W2mDBpcqLgRM7PKQqq2tparFy5ElarFc899xxSUlIAAJWVlXjppZewefNmjB8/HgkJCZ26Xn5+PlavXg2VSoXly5cjPj5e2v/iiy9i9erVGDduHMLCHCs/BwYGYvr06Zg2bRq0Wq20/+TJk3jjjTewefNmjB07FmPGjOnFd+7egiNCMcp4Cec09lW19527zYKKiIiIyIntvV4JqOwF1URrMTy8WVDdiVOMv9qzZw/q6+uRnJwsFVMA4OvriyVLlgAANm3a1OnrbdmyBRaLBXPmzJGKKQAIDw/H/fffD4vFgq1bt7Y455lnnsGcOXNaFFMAkJSUhJkzZwIADh061OX3NtClhaul7YMWfxgNjQLTEBEREVF7asorcULRbO2pOD+BaVyHUxRUmZmZAIDU1NRWryUlJUGlUiErKwtGo7FT1zt58mS715s0aVKLn9kZMTH2efcrKio6fQ7ZTZ4yGmqLCQBQq9Qj88gZwYmIiIiIqC0HD2fBLLcPYPM21WFsymjBiVyDUxRUubn2KbVjY2NbvaZUKhEdHQ2TyYT8/PxWr39dXV0dSktLAdgnofi6gIAAeHl5oaSkBPX19Z3KV1RUBADw8fHp1PHk4OHliRQ4FqLak1PTwdFEREREJMreIqu0PV1TBZXKKZ4OcnrC/5Tq6+tRV1cHAPD392/zGH9/f1y7dg2lpaVtFknNNRVTHh4erYbvNQkICEBNTQ1KS0sRHR3d4fXq6uqQnp4OAJgwYUKHxzb37LPPttqnVqvx2muvAbA/syWSUmn/qw8K6vtpMBeNi8GBs/ZeqkxlKORmKwLCQu5wFonUn/cHuRbeG9Qe3hvUHt4bruHGhcvI1jo+n907fXS//J25w/0hvIfKYDBI2xpN29MyNu1vfuydrtfetbp6vX/961+orq7G0KFDMXHixDseT62lTp8AP1MtAMAiV2DbzmOCExERERFRc5vTz0rbMcZyjBzDicQ6q1d6qN58803cunWrS+c89dRTiIuL640f34LNZgNgXwepp7788kscPnwYnp6eePrpp7t0zRUrVnT4emlpqZRVhKZvAUpKSu5wZO+YrqvBerMnAGBbvhFziot75e+I+kZ/3x/kOnhvUHt4b1B7eG84P4vFjJ3VWuCrucRm+FukUV99TfT9IZPJWsz83R29UlCVlJR06vmm5hob7bO9NR+W19jYCL1e36lj26PT6QB03PvUmevt27cPH3/8MTQaDX75y18iJIRD1Hpi5vghWH/M/szaNW0Ibl26hugRvV9QExEREVHXXDh+DiVq+1wBcpsVaZNHCU7kWnqloGp6Lqg79Ho99Ho96uvrUV5e3mZBVV5eDqBzzx01HVNXVweDwdBm0VRWVtbh9Y4fP46///3vUCgU+NnPftZi6nXqnsFx0Rh84CBuqO1/5ntP3cBjLKiIiIiIhNubXQwo7fMKjDYXIyBopOBErkX4M1SAY1ry69evt3rNbDYjNzcXKpUK4eHhd7yWh4eHVCjl5OS0er2srAw1NTUIDAxss3g7f/483n77bQDA008/zYV8e9HMIMfttq/RB+ZOToNPRERERH3DUF2DwwiW2nfFeApM45qcoqBKSkoCABw9erTVaydPnoTJZEJCQgLUanWr17t6vSNHjrQ4prnr16/jj3/8I8xmM/7nf/6nzXWsqPumTxkFuc0CAChXe+Ncxtk7nEFEREREfenoodNoUNpHdOksjUiZlCg4ketxioJq1qxZ0Ol0OHHiBI4dc8wAV1VVhQ8++AAAsGjRolbnLVu2DMuWLZOGBDZZsGAB5HI5du7cicuXL0v7CwoKsG7dOsjlcixYsKDFOfn5+fjDH/6AhoYGPP7440hLS+vFd0gA4OfngyRLszWprpR3cDQRERER9bW9+Y4RQ5MVFdBq258pm9omfB0qAPD09MSTTz6Jt956CytWrMDIkSPh5eWFrKws1NXV4e6770ZiYutquWkiDLPZ3GJ/eHg4lixZglWrVmH58uUYPXo0FAoFzp49C6PRiKVLl7YaPvj222+juroa3t7euH79Ov7617+2+nkRERG47777eu+ND0AzY31xwr6OM44qQlFfUQm9n6/QTEREREQDUWnubZzVOGa4uysxQmAa1+UUBRUApKam4uWXX8batWtx5coVmM1mREREYN68eZg5c2aXr7do0SKEhoZi48aNuHjxIgAgNjYWixcvRnJycqvjmxYXrq6uxv79+9u85siRI1lQ9dCElFHwuH4OdUodGhVqHD50BrMXzRAdi4iIiGjA2X/0IqyyKABAsKkaIxJaf0amO3OaggoAhg8fjhdeeKHTx69Zs6bD15OTk9ssntrSVo8U9T6NWoUp6irssNqnt99baMZswZmIiIiIBhqrxYK9NVrgqwmx07wboZA7xdNALod/atTv7hoTLW2f00Wg6FqOuDBEREREA9D1MxdwSxsktWemDBOYxrWxoKJ+N3zEYISZqqT2/uOXOziaiIiIiHrb3qw8aXuYqRThEcEdHE0dYUFF/U4mkyHNzzGRyN46D1i/NrEIEREREfUNU3090m2BUvuuCM7s1xMsqEiItNQR0na+NgBXMrMEpiEiIiIaOE4ePoVqlX0BX5XVjCmTEwQncm0sqEiI0BB/jDI3W5PqYpHANEREREQDx56cWml7oqwMXh46gWlcHwsqEmZmtIe0fUAWgsaaGoFpiIiIiNxfVUEhTmgc603NHB4iMI17YEFFwkxJHQWNxb46d51Sh4xDpwQnIiIiInJv6YfPwyy3r5zka6rDuHHxghO5PhZUJIxep8EkZYXU3n3bKDANERERkXuz2WzYU+FYhjbNsw5KBcuBnuKfIAk1a1S4tH1GG47SnFsC0xARERG5r5ysi7iucwzxmzVhqMA07oMFFQk1KnEIgk3VAACrTI69GdmCExERERG5p91nHF9cDzWVIjomTGAa98GCioRSyOW4y6dRau+t0cFqsQhMREREROR+TAYD9luarT0VphKYxr2woCLhZjZbkypPG4Dsk+cFpiEiIiJyPycOn0a1yj7DsspqwjSuPdVrWFCRcKFhgUgwOdah2nOhUGAaIiIiIvfTfO2pFFspvLw8OjiauoIFFTmFu6L00vZBBMNQW9vB0URERETUWRVFJchUOZ6XmjU8SGAa98OCipzC5MmJ0Frsz1LVK7U4euis4ERERERE7mH/4fOwyBUAAH9TDUYnjbjDGdQVLKjIKeh0WkyRl0nt3XmNHRxNRERERJ1htVpbrD01U18HpVIhMJH7YUFFTqP5mlRZ6hAU5+YLTENERETk+q6fv4qbmmaz+02IE5jGPbGgIqcxYkw8woyVAAAb16QiIiIi6rE9Z3Ol7WGNRYgcHCkwjXtiQUVOQy6XY6a3QWrvqdbBarUKTERERETkuowGI/ab/aX2XWHKDo6m7mJBRU5lZuoIyGz2IqpQ44sLpy8JTkRERETkmo4fPYNapX0mZbXFiKlTxghO5J5YUJFTCY4IwWhjszWpzvM5KiIiIqLu2H3DsQxNqq0Ent6eAtO4LxZU5HTuitJJ24dsQWjgmlREREREXVJWWIxTqmCpPWt4YAdHU0+woCKnkzp1DPRm+7NUBoUGhw6cEZyIiIiIyLXsO3weVpl9evRAYzUSkxMEJ3JfLKjI6Wh1OkxVlkvt3flGgWmIiIiIXIvVYsHuSo3UnulZD4WCa0/1FRZU5JRmj42Wti9ow5B3JUdcGCIiIiIXcunUReRpHLP7zZo0XGAa98eCipxS/MghiDaWSe1dx68KTENERETkOnadL5C2E41FCIsMFZjG/bGgIqckk8kwO8gmtfc2+sDcyKF/RERERB2pr6zGIZljMorZMR4C0wwMLKjIaaVNSYTSagYAVKi9cPLIKcGJiIiIiJzbwYOnYVDYn5/Smw1InZwoOJH7Y0FFTsvHzwcTbCVSe+f1aoFpiIiIiJzfriKrtD1dXQGtVtPB0dQbWFCRU5sz3NFlnakKR0V+ocA0RERERM4r9+JVZGsdz0vNGT9YYJqBgwUVObUx40cgwFQDALDIFdh7+JzgRERERETOadfJ69L2YGMZ4uIHiQszgLCgIqemVMhxl1eD1N5VpYPVYhaYiIiIiMj5mBoM2GdqNlV6KD/m9xf+SZPTmzVphLSdpw1A9gn2UhERERE1d+LQSVSpPAEAKqsJM6ZwMor+woKKnF5YeBASTcVSe+eFIoFpiIiIiJzPrpt10nYKyuDt7SkwzcDCgopcwuzBjl8KhxShqK+oEJiGiIiIyHmU3byFk5oIqT1nFBfy7U8sqMglpKYmQm82AAAMCg0Opp8WG4iIiIjISew5ehFWmf1jfZCpGolj4wUnGlhYUJFL0GpUmK6pktq7Smyw2WwCExERERGJZzWZsLvWMZJnlq8RCjk/4vcn/mmTy5iTHCttZ+vCcet8tsA0REREROJdyDiNAq19dj+ZzYpZk0YKTjTwsKAilxEXF4XBpnKpvetUjrgwRERERE5gZ3aZtD3GUorgEP8Ojqa+wIKKXMqsMKW0vc8cCFN9XQdHExEREbmvupISHFaGS+3Zcb7iwgxgLKjIpcyYkgCV1b6wb5XaExkHTgpORERERCTG/oNnYVSoAQCe5gakTORwPxFYUJFL8fbUI1Xu6NreecsgMA0RERGRGDarFTvLHCN30vS1UKuUHZxBfYUFFbmcuYmOru3T2ggUXcsRF4aIiIhIgGunz+O6LkRqz50YJzDNwMaCilxOYmIcwkz2KdRtMjl2HbssOBERERFR/9qRlS9tDzOVIiYmTGCagY0FFbkcmUyG2f5mqb270QfmRqPARERERET9p6GqGulo1jsVpRWYhlhQkUuaNTUBCqsFAFCm9sGpI6cEJyIiIiLqHwcPnkaD0l5E6S0GTJmSKDjRwMaCilySn78PJthKpPaO6zUC0xARERH1n52FVml7mqoSOq1GYBpiQUUua+7wIGn7hDoMZXmFAtMQERER9b2bF64iWxsqteeOHyQuDAFgQUUubMz4EQg0VgMArDIFdh8+LzgRERERUd/aceqGtB1rLEVc/CBxYQgACypyYUqFHLN9GqT2rhoPWCzmDs4gIiIicl2N9Q3YZwqQ2nNCFQLTUBMWVOTSZk8aCbnNPo64SOOLrIxzghMRERER9Y0jh06jVqUHAGgsRsyYOlpwIgJYUJGLCwoLwjhTkdTekV0mMA0RERFR39l5q1HaniIrhYeXh8A01IQFFbm8uUO8pe1jihBUlpYLTENERETU+25fv4VzmmaTUYwOF5iGmmNBRS4vedIY+Jns06ab5UrsPchhf0RERORedmZckbajG8swLDFeYBpqjgUVuTylSom7dI51qHaWq2C1Wjs4g4iIiMh1GBtN2GPwkdpzAi2Qy/kx3lnwb4LcwpwUx7c0eRo/XDh1QWAaIiIiot5z/PBpVKvsz0uprCakcTIKp8KCitxC2KBIjG4skNo7z3GRXyIiInIPO244RuJMshbD299XXBhqhQUVuY25sV7S9mFZCKo5OQURERG5uMIbt3BG3WwyisQwgWmoLSyoyG2kTBoDb1MdAMCoUGHvgTOCExERERH1zI6jl2CT2T+yhxsrMGrscMGJ6OtYUJHbUGtUmKV3dIlvL1fDarEITERERETUfUaDAbsa/aT2XE5G4ZT4N0JuZe4kx7c2edoAXDieJTANERERUfcdO3ASVSpPAPbJKO6aNkZwImoLCypyK+FRoRhjKpLa2y4WC0xDRERE1H3bcw3S9mSUwsfXq4OjSRQWVOR25sc51mk4ogxHZQGLKiIiInItt7OvIUsbLrXnj40UmIY6ohQdoLns7GysXbsWly9fhtlsRmRkJObNm4e0tLRuXS8zMxMbNmxATk4OAGDQoEFYvHgxxo8f36nz9+/fj7/+9a8AgEcffRT33Xdft3JQ/5qQmgjfq6dQqfKEWa7E7kNn8eBDs0XHIiIiIuq0HcevAbJoAECUsRwjEiYJTkTtcZoeqoyMDCxfvhynT59GTEwMxo4di8LCQqxcuRLvv/9+l6+3ZcsWvP7667h8+TKGDRuGUaNG4dq1a3j99dexZcuWO55fXV2NVatWQSaTdeftkEAqpQJzvOql9s5qPSwmk8BERERERJ3XWFeHPeYAqT0vVMbPpE7MKXqoamtrsXLlSlitVjz33HNISUkBAFRWVuKll17C5s2bMX78eCQkJHTqevn5+Vi9ejVUKhWWL1+O+Ph4af+LL76I1atXY9y4cQgLa38e//fffx8GgwFTp07FgQMHev4mqV/NmTISn+8qhk0mR4HGH1lHT2PstAmiYxERERHd0eH0k6hRBQEA1BYT0qZyMgpn5hQ9VHv27EF9fT2Sk5OlYgoAfH19sWTJEgDApk2bOn29LVu2wGKxYM6cOVIxBQDh4eG4//77YbFYsHXr1nbPP3v2LA4cOIAHH3wQISEh3XhHJFpIaCCSLCVSe9vlCoFpiIiIiDrHZrNhR55jZM1UZRm8vPQCE9GdOEVBlZmZCQBITU1t9VpSUhJUKhWysrJgNBo7db2TJ0+2e71Jkya1+JlfZzQa8a9//QsRERFYvHhxp34eOad5wwOl7Qx1BMpv5QlMQ0RERHRnueeycUHXbDKKpBiBaagznKKgys3NBQDExsa2ek2pVCI6Ohomkwn5+fl3vFZdXR1KS0sB2Ceh+LqAgAB4eXmhpKQE9fX1rV5fs2YNioqK8IMf/ABKpVOMiKRuGp88AoEm+0K/FrkCuw6dF5yIiIiIqGPbT+ZI24NN5YgfNkhYFuoc4RVDfX096urqAAD+/v5tHuPv749r166htLS0zSKpuaZiysPDA1qtts1jAgICUFNTg9LSUkRHR0v7c3JysHnzZqSlpWHkyJHdeDcOzz77bKt9arUar732GgAgMDCw1ev9qalYDAoKEpqjry0ItmHVV6P9dtR740fe3lBqNGJDuYCBcn9Q1/HeoPbw3qD28N7ovPqKCuyF43GTe4d4ITg4WGCivucO94fwHiqDwbFgmaadD7pN+5sfe6frtXet9q5ntVrxj3/8A3q9Ht/5znfuHJxcwv3zUyG3WQAAJRpfHNyeLjgRERERUdu2bUlHvVIHANBaGrFgHqdKdwW90kP15ptv4tatW10656mnnkJcXFxv/PgWbDYbAHR5asktW7bg2rVrePLJJ+Hl1fNVqFesWNHh66WlpVJWEZq+BSgpKbnDka5NLgcm2EpxTGb/tmfd+WKMmuTe77k3DJT7g7qO9wa1h/cGtYf3RufYbDasz20E7PUUpqsr0VBfi4b6WrHB+pjo+0Mmk3U483dn9EpBVVJS0qnnm5prbGwEgBbD8hobG6HXt57FpK1j26PT2e/Cjnqzvn69kpISfPrppxgxYkS3FxEm5zV/VCiOXbQXr5naSJRcy0HQkEFiQxERERE1c/3UOVzROT7Yz58wRGAa6opeKaiangvqDr1eD71ej/r6epSXl7dZUJWXlwPo3HNHTcfU1dXBYDC0WYSVlZW1OPb8+fNobGxEdXU1Xn755RbHNlXLu3btwunTpzF8+HA88sgjXXiHJNqYsUMRcvYEilTesMrk2Hn0Eh5lQUVEREROZPuZ24B6MABgqKkMQ4YMF5yIOkv4pBQAEBMTg4sXL+L69euIjIxs8ZrZbEZubi5UKhXCw8PbuYKDh4cHAgMDUVpaipycHAwf3vJmLCsrQ01NDQIDA1sVb3l5ecjLa3tq7eLiYhQXF7dZ8JFzU8jlmBNowQdV9vZOoz++0dAA1Ve9mUREREQi1ZWWYb/c8Tl3Xgw/o7gSpyiokpKScPHiRRw9ehTTp09v8drJkydhMpkwbtw4qNXqTl9vx44dOHr0aKuC6siRI9IxTdLS0tod6rdmzRp8/vnnePTRR3Hfffd1/k2RU5k9NRGfbMyBWa5EudobGeknMGXeNNGxiIiIiLAv/TQMyggAgN5iwLTJiYITUVcIn+UPAGbNmgWdTocTJ07g2LFj0v6qqip88MEHAIBFixa1Om/ZsmVYtmyZNCSwyYIFCyCXy7Fz505cvnxZ2l9QUIB169ZBLpdjwYIFffRuyBn5+XpisrxMam+91blFoomIiIj6ktVixtYKR6fBLF01tBqVwETUVU7RQ+Xp6Yknn3wSb731FlasWIGRI0fCy8sLWVlZqKurw913343ExNaVetNEGGazucX+8PBwLFmyBKtWrcLy5csxevRoKBQKnD17FkajEUuXLu3U8EFyL3ePjUL6aXshlaWLQO6Fy4geGS84FREREQ1kFzLO4pbWsQbT/NRhAtNQdzhFQQUAqampePnll7F27VpcuXIFZrMZERERmDdvHmbOnNnl6y1atAihoaHYuHEjLl68CACIjY3F4sWLkZyc3NvxyQWMGDkYMRlHcFNtX0B6W2YOfsiCioiIiATacrEU0EQDAEabihAZxckoXI3MJnIxpAGsoKCA61AJsHX7Mfy91AcAoDcb8J/7Y6Hz9hacyvkM1PuD7oz3BrWH9wa1h/dG+8rzC/HE7jJY5AoAwC8HN2LS5DGCU/Uv0fdHb6xD5RTPUBH1lxnTxkBnsa9DVq/UYv/+02IDERER0YC18+B5qZjyN9VgwsQEwYmoO1hQ0YCi12sxU+WYxGRbsQxWq1VgIiIiIhqIzEYTttd5Se15Pg1QKhUCE1F3saCiAefuiXHS9g1tELJPXxCYhoiIiAai44dPoUxtf+xAYbVgzlT2TrkqFlQ04EQPiUKCsUhqb8kqFJiGiIiIBqKt12ul7VRrEQKC/AWmoZ5gQUUD0oLBemn7sDwUFSVlHRxNRERE1HtuX8/FGY1jCZ8FCaEC01BPsaCiAWni5LHwM9UAAMxyJXYdPCc4EREREQ0U245dkbajGsswchynSndlLKhoQFKpVZird3S1b6/Uwmy2CExEREREA4GhvgF7jI7hfXcHWyGX8yO5K+PfHg1Yc6eOgNxmL6JK1D44eeys4ERERETk7tLTT6FOqQMAaC2NSJs2Vmwg6jEWVDRgBYaHYqLJMSHF1uwKgWmIiIjI3VmtVmzNd4yISZOXwMPHq4MzyBWwoKIBbcGoYGn7lDoU+ddzBaYhIiIid3blzCVc1wRJ7fkThghMQ72FBRUNaInJoxBhtPdM2WRybDuSLTgRERERuastZ/Kk7ZGNRRg8bLDANNRbWFDRgCaXyzE/2Cq1d5kDYaip7eAMIiIioq6rKCjGQUWY1L471lNgGupNLKhowJs1Yxy0lkYAQJ1Sh337MgUnIiIiInez48BZmOVKAICfqRaTpo4RnIh6CwsqGvA8PPWYqXJMSLGlWA6rhVOoExERUe8wNTZiW5231J7nUw+VUikwEfUmFlREABakDpW2b2qDcD6DU6gTERFR7zi6/wTK1faCSmk1Y9700YITUW9iQUUEIHpwBEabiqT25oulAtMQERGRu7DZbNiS2yi1J6ME/gG+4gJRr2NBRfSVRcMcq5YfU4ejOOeWwDRERETkDnKyLuKCLlxqLxwfIzAN9QUWVERfGT9xFIKNVQAAq0yBbYcuCk5ERERErm7TSccXtEOMZRg2glOluxsWVERfUSrkmB/omIxipykAjbV1AhMRERGRK6suKkG63DFV+sIYDWQymcBE1BdYUBE1MydtDNQWEwCgWuWBg/tOCE5ERERErmpX+mkYFWoAgLe5HlM5VbpbYkFF1Iy3lwemK8uk9uZCcAp1IiIi6jJzYyO21nhJ7TmeddCoVQITUV9hQUX0NQubTaF+TReCKyc4hToRERF1TeaBEyjW+AIA5DYL5k9PEBuI+gwLKqKviY2NwEhTidTedL5YYBoiIiJyNTabDZtzGqR2iq0UwUF+AhNRX2JBRdSGhUN9pO3DqgiU594WmIaIiIhcye3z2Tiji5TaC8dGdnA0uToWVERtSElNgL+pFgBgliux/eB5wYmIiIjIVWw+cUPajjGVIyEhVmAa6mssqIjaoFLIMd/fKLW3N/rBVM8p1ImIiKhjdSUl2CtzLOS7IIpTpbs7FlRE7Zg7fTSUVjMAoELtjcOcQp2IiIjuYPe+UzAoNQAAD7MBM6aOFpyI+hoLKqJ2+Pl6YqqiXGpvzLfCZrUKTERERETOzNxoxOZqT6k927MWOg2nSnd3LKiIOnDPRMeY5yu6MGRzCnUiIiJqR+aBEyjU+gMA5DYrFk4dKTgR9QcWVEQdiIuLbDGF+sZznEKdiIiIWrPZbNjYYqr0EoSE+AtMRP2FBRXRHSyK95a2D6sjUMIp1ImIiOhrcrKykaWLkNqLxkZ0cDS5ExZURHeQkpKIIFM1AMAqU2DrwYuCExEREZGz2XQqV9qONZZiVMIQgWmoP7GgIroDpUKOBQEmqb3dGABDLadQJyIiIrvKolLsl4dJ7XtiOFX6QMKCiqgT5swYA63Fvi5VrUqPfftOCk5EREREzmJ7+lmY5PbZ/HxNtZg6dazYQNSvWFARdYKXtyfSlKVSe1OxHFZOoU5ERDTgGRuN2FrrJbXne9VBreZU6QMJCyqiTrpn0lBp+5YmAGczsgSmISIiImdwOP0kKtT2gkppNWP+dC7kO9CwoCLqpMjBUUgy5kvtDRfLOziaiIiI3J3NZsPGW0apPR1F8AvyE5iIRGBBRdQF9wx3rCeRqQ5D3g1OoU5ERDRQXcq6jKuaYKm9aHyMwDQkCgsqoi4YOzERkY1lUnvTkcsC0xAREZFIG085vlgd1ViIISPjBKYhUVhQEXWBXKHAomDHZBR7zAGoqagWmIiIiIhEKL5dgCOKUKl9T6yHwDQkEgsqoi6aOTMZnuYGAIBBocGuvZmCExEREVF/25J+HlaZAgAQbKzChKnjBCciUVhQEXWR1kOHObpKqb2lUgez0dj+CURERORWGqprsNPkeK56QYAJSqVSYCISiQUVUTcsmJ4Iuc0CACjW+OLo/hOCExEREVF/2bMnE7VKPQBAazFidhp7pwYyFlRE3RAcGojJtmKpveGmETabTWAiIiIi6g9mowkbyzVSe5a6HF7efH5qIGNBRdRN9yUPkrazdaG4dOKsuDBERETUL44fOIECjX2tKbnNisXTRwpORKKxoCLqpqEjBmOkydFLtT6rSGAaIiIi6ms2mw0bbjRI7RRrMULDgzs4gwYCFlREPXDvMMdq6MfUkSi8ekNgGiIiIupLV06ewwVduNS+d3yUwDTkLFhQEfXAhJRRCDNWAgCsMjk2HuZCv0RERO5q/ZkCaXuYqQQjRg0RmIacBQsqoh5QyOW4J1wmtXfZQlFTUiowEREREfWF4hs3cVgdIbXvHeotMA05ExZURD1014xxjoV+lRrs2HtKcCIiIiLqbZsOXpIW8g0xVSElNVFwInIWLKiIekinVWOeV63U3lzrDVNDQwdnEBERkSupKyvDTmuI1F4YYoNSwY/RZMc7gagXLEwbDaXVDAAo0/jg0J4MwYmIiIiot+zccxL1Si0AQG82YA4X8qVmWFAR9YIAfx9MU5RJ7fUFgNViFpiIiIiIeoO5wYBNNV5Se65HDfQ6TQdn0EDDgoqol9w7KU7avq4LwfnDJwWmISIiot5wZG8GSjS+AACF1YJFaXx2ilpiQUXUSwYPjsAYc4nUXp9dITANERER9ZTVYsH6fKvUniIvQ1Cgr7hA5JRYUBH1onsTHauln9BG4fa5SwLTEBERUU9cOnYSV3ShUntxaqzANOSsWFAR9aKkcfGIMtl7pmwyOTYcvyE4EREREXXXlxfKpe1R5lIMHRIpMA05KxZURL1IJpNhcYzjQdW9ikhU3M4XmIiIiIi64/aFbGRoo6T2vSMDBaYhZ8aCiqiXzZg6Gr7mOgCAUaHClv1ZghMRERFRV63PuAGbzP5ROcJUiQnJwwQnImfFgoqol2lUSizyb5TaW0xBaKiqEpiIiIiIuqIirwB7FRFS+75oJeQymcBE5MxYUBH1gflpY6G12IuqWpUeu3ZnCk5EREREnbVpfxZMchUAwNdUi7RpYwQnImfGgoqoD3h56TFHWym1N1TqYTYaxQUiIiKiTqmvrsY2k+N5qUX+jVCrVAITkbNjQUXURxZPHwWF1QIAKNb44tC+E4ITERER0Z3s3J2JWqUeAKC1NGL+zLFiA5HTY0FF1EeCQwMxDcVSe90tC6xWawdnEBERkUgmoxEbKvVSe56mAl5eHgITkStQig7QXHZ2NtauXYvLly/DbDYjMjIS8+bNQ1paWreul5mZiQ0bNiAnJwcAMGjQICxevBjjx4/v8LysrCxs3boVV65cQV1dHby8vDBo0CDMmTMHycnJ3cpCA9N9qbHYl9EAALihDcLZjCyMTeU4bCIiImd0cH8mStV+AACF1YJF00cJTkSuwGkKqoyMDKxYsQI2mw0jRoyAl5cXzp07h5UrV+LmzZt47LHHunS9LVu24L333oNCoUBiYiKUSiXOnj2L119/HY8//jgWLFjQ5nkffvgh1q9fD6VSiWHDhsHHxwcVFRW4cOEC/Pz8WFBRlwweGoNxB/filDoMALDuQjnGpgoORURERK1YrVZ8ecsCfLWc5DQUIjiMBRXdmVMUVLW1tVi5ciWsViuee+45pKSkAAAqKyvx0ksvYfPmzRg/fjwSEhI6db38/HysXr0aKpUKy5cvR3x8vLT/xRdfxOrVqzFu3DiEhYW1OG/Hjh1Yv349hgwZgueeew6BgY4HEhsbG1FUVNRL75gGkvtHBeDUFfv2aU0Yrl26jiHDY8WGIiIiohZOHz+HHI3js999E/n/auocp3iGas+ePaivr0dycrJUTAGAr68vlixZAgDYtGlTp6+3ZcsWWCwWzJkzRyqmACA8PBz3338/LBYLtm7d2uKcuro6fPjhh9DpdPjFL37RopgCAI1Gg+jo6O68PRrgEpMTMKTR8SzVlxk54sIQERFRm9ZdKJe2k4z5GDxssMA05EqcoqDKzLSv0ZOa2nosVFJSElQqFbKysmDs5LTTJ0+ebPd6kyZNavEzmxw6dAgNDQ2YMmUK/Pz8upSfqCNyuRz3R6ul9kF5CIpu5QtMRERERM1dPX8VZ9WhUvv+kYEdHE3UklMUVLm5uQCA2NjWXatKpRLR0dEwmUzIz7/zh9C6ujqUlpYCsE9C8XUBAQHw8vJCSUkJ6uvrpf1ZWVkAgNGjR6OyshKbNm3CP//5T6xevRoZGRmcnY16ZNKMZIQYKwEAVpkCG9IviA1EREREknXHc6TtIY3FSJjQucdMiAAneIaqvr4edXV1AAB/f/82j/H398e1a9dQWlraZpHUXFMx5eHhAa1W2+YxAQEBqKmpQWlpqTSM7/bt2wCAkpIS/P3vf29RbG3cuBGDBw/G888/327Gr3v22Wdb7VOr1XjttdcAoNWQwv6mVNr/6oOCgoTmGEgejlLgz189hrfTEownrTb4hQSLDdUO3h/UHt4b1B7eG9QeZ783ci9fw2FluNT+dkIAQkJCBCYaWJz9/ugM4T1UBoNB2tZoNG0e07S/+bF3ul5712rverW1tQCAjz76CKGhofj973+P999/H7/73e8wePBg3LhxA2+++SZsNtsdMxC15b7FafAy2wv1RoUan36ZLjgRERERfbj1OKwy+0fiUFMVZs+bJjgRuZpe6aF68803cevWrS6d89RTTyEuLq43fnwLTQWPTCbr0nlNQ/rUajVeeOEFeHt7AwDi4+Pxwgsv4KmnnsKVK1eQlZWF0aNH3/F6K1as6PD10tJSocVZ07cAJSUlwjIMRAu9a/FJvX3BwC+rPbHw2jVov7rXnAnvD2oP7w1qD+8Nao8z3xsV+UXYbg4GFPb24lAbKirKOz6JepXo+0Mmk7Wa+bureqWgKikp6dTzTc01NjYCQItheY2NjdDr9Z06tj06nQ5Ax71ZbV1Pp9OhpqYG48ePl4qpJj4+PkhKSsKRI0dw4cKFThVURG1ZcNc4fLn+BgwKNapVHtix6zgWPzBLdCwiIqIBaeO+MzAq7I9/+JjqMCstSXAickW9UlA1PRfUHXq9Hnq9HvX19SgvL2+zoCovt39T0JnnjpqOqaurg8FgaLMIKysra3W9oKAgFBcXtzt+s2l/VVXVHTMQtcfHxwtztBXYaLKPzV5f6Yn5DQaodXf+soCIiIh6T11ZObaagqRPw/f4G6DVqjs+iagNwp+hAoCYmBgAwPXr11u9ZjabkZubC5VKhfDw8Favf52Hh4dUKOXk5LR6vaysDDU1NQgMDGxRvA0ebF9roOlZqq9r2t+ZXjKijtw7czSUVjMAoFTjg/RdRwUnIiIiGni27s5EvdI+sklvMeDuWeydou5xioIqKcl+Ax892vqD5cmTJ2EymZCQkAC1unPfGnR0vSNHjrQ4pklycjIA4MKFC62mSLdarbh48SIAR+FF1F1BQX6YoSiT2muLVbCYOrfGGhEREfWcoaYGG+sc647e7VEDTw+dwETkypyioJo1axZ0Oh1OnDiBY8eOSfurqqrwwQcfAAAWLVrU6rxly5Zh2bJl0pDAJgsWLIBcLsfOnTtx+fJlaX9BQQHWrVsHuVyOBQsWtDhn5MiRiI+PR15eHtauXdvitc8++wwFBQXw8fHBxIkTe/x+iR6YNhwym71wz9MG4OjeDMGJiIiIBo49u46jUu0JAFBbTbjnrrFiA5FLE74OFQB4enriySefxFtvvYUVK1Zg5MiR8PLyQlZWFurq6nD33XcjMTGx1XlNE2GYzeYW+8PDw7FkyRKsWrUKy5cvx+jRo6FQKHD27FkYjUYsXbq0zeGDTz31FF588UWsWbMGhw4dQmRkJG7fvo28vDyo1Wr85Cc/4ZA/6hWRUSFIxSUcgf1Zqi9yLZhkMUOucIp/kkRERG7LbDDgywod8NUKO7PUFfDza/05k6iznObTW2pqKl5++WWsXbsWV65cgdlsRkREBObNm4eZM2d2+XqLFi1CaGgoNm7cKA3Xi42NxeLFi6XhfV8XGhqKN954A2vWrMGpU6dw4sQJeHp6YsqUKXjggQcQFRXVo/dI1NxDKbE4kmFf1PqaLgRnD2Zi7IwUwamIiIjc24Hdx1CksX+hKbdZcF/aKMGJyNU5TUEFAMOHD8cLL7zQ6ePXrFnT4evJycntFk/t8fX1xQ9/+MMunUPUHXFDozD2cDpOK4MBAJ9fqcGY6bYur6FGREREnWMxmbC2QA589bjUdHkpQkNYUFHPOMUzVEQD1YNJEdJ2li4Sl4+fFheGiIjIzZ3Yn4FcnWOJnAemxAtMQ+6CBRWRQIkJsYg3lUrtz88WC0xDRETkvqwWC7640Si1J1qLEBMTJjARuQsWVEQCyWQyPDQqQGpn6GKQe/aCwERERETu6cKRTGTrHZOSPTRxkLgw5FZYUBEJNiF5OKJMFVL7i+O5AtMQERG5H5vNhs8vVkrtBHMJhg2LEReI3AoLKiLB5DIZHhyil9rpmigUXr4mMBEREZF7uZ55Bqf00VL7obGhAtOQu2FBReQEpk1ORLCpGgBglSmw7vDlO5xBREREnbXmdJG0PcRUjrGj4wSmIXfDgorICSgVcjwQ6fjnuEsRhbKbtwQmIiIicg83sy7hqM4xvO8bI3y4RAn1KhZURE7irhlj4G+qBQCY5UqsS+fkFERERD31WbNnk6NNFUiZOFJgGnJHLKiInIRGpcJ9IWapvd0WhooCTqNORETUXXmXr+GQOlJqf2OIDnL2TlEvY0FF5ETmzRwHH1MdAMCoUGP93jOCExEREbmuz49cg1Vm/7gbbqzE5MmjBScid8SCisiJaLUa3OvXILW3moNRXVbRwRlERETUlsKbedgvj5DaD0bJoVTwoy/1Pt5VRE7m7lnj4WmuBwAYFBps3H1KcCIiIiLXs/bARVjkCgBAsLEKM9KSBCcid8WCisjJ6D11uMezWmpvbvBDbXWtwERERESupbSwGLttjrWmHgixQKVUCkxE7owFFZETWjgrCXqzAQBQp9Rh6+5MwYmIiIhcx7q952CW2wsof2MNZs1KFpyI3BkLKiIn5OXrjQXaMqm9ocoTDXX1AhMRERG5hoqScuwwB0nt+wPqodaoBSYid8eCishJLb5rHLSWRgBAtcoD23YeF5yIiIjI+a3ffQpGuQoA4GOqxdzZEwQnInfHgorISfkE+WOeqlRqf1npAUM9e6mIiIjaU1Vajq3GQKl9r08ttHq9wEQ0ELCgInJi980aC7XFBACoVHli145jghMRERE5r407T8Kg0AAAPM31uHvuRMGJaCBgQUXkxPyDAzBH7eilWleuh5HPUhEREbVSU1aGzcYAqb3YuxZ6D/ZOUd9jQUXk5O6bORpKqxkAUKrxwa4dRwUnIiIicj6bdmSiXqkDAOgtBiyYPV5wIhooWFARObngkADMVjp6qb6o8ICxrk5gIiIiIudSU1qGDUbHzH73eNXAy8tDYCIaSFhQEbmAB2eNYS8VERFROzbtPOHonTIbcA97p6gfsaAicgHBwX6YrXKsS/V5pSd7qYiIiADUlJRigzFYat/jXQsvLz47Rf2HBRWRi3jwLsezVGVqH+zafkRwIiIiIvE27sxs2Ts1J0lwIhpoWFARuYjWvVReMNayl4qIiAaumpISbDQ1653yqYOXJ3unqH+xoCJyIQ/OatZLpWEvFRERDWwbd5z82rNT4wQnooGIBRWRCwkO8sNsdbnU/rzKC8aaWoGJiIiIxKgpLsFGs6N3arEve6dIDBZURC7mwbsSW/RS7dzBXioiIhp4Wjw7ZTFg0Ww+O0VisKAicjHBQX6YrWnWS1XtA2NNjcBERERE/aumqBgbzSFSe7FPHbw8dAIT0UDGgorIBT3UbMa/crU3dm7nulRERDRwbNjleHbKw8zeKRKLBRWRCwoK9MWc5r1UNT5oZC8VERENADVFRdhkDpXa9/jVs3eKhGJBReSims/4V672xs5t7KUiIiL3t2HnSdQrtQC+6p2axd4pEosFFZGLCgrwxRxNhdT+otYHhmr2UhERkfuqKijEJkuY1L7HrwFeHlqBiYhYUBG5tIdmjYaqWS/VNj5LRUREbuzL3adb9E7dM4vrTpF4LKiIXFhggA/max3PUn1R54f6qiqBiYiIiPpGeX4hNlvDpfb9/g3wZO8UOQEWVEQu7sFZY6GxGAEA1SpPbNp+XHAiIiKi3vfF7rNoVKgBAN7meizkzH7kJFhQEbk4P39vLNQ7nqX6siEANRWV4gIRERH1suJb+dgGR+/Ug4EG6HUagYmIHFhQEbmB++ckQW82AADqlDqs35EpOBEREVHv+WzfBZjlSgCAv6kGd89OFpyIyIEFFZEb8PbxwmLPSqm9sTEQVaUV7Z9ARETkIgpu5mE3HDP7fSPYBI1GLTARUUssqIjcxD1zk+FprgcAGBQarN11SnAiIiKinvs0/RIscgUAINhYhdmzJghORNQSCyoiN+Hp5Yn7vRwz/G0xBqKsuExgIiIiop65deMW9sscvVPfDDNDrVEJTETUGgsqIjeycO4E+JpqAQBGhRpf7DotNhAREVEPfJJ+GVaZ/eNqeGM50u6aKDgRUWssqIjciM7TEw/61knt7ZYQFOfmC0xERETUPdfPX8VBZYTUfiRKDqWavVPkfFhQEbmZefNTEGCsBgCY5Up8uvec4ERERERd99GxHGk7prEMU+/is1PknFhQEbkZjVaLh8MsUnuPIgL5V64LTERERNQ12SfO4rgmUmo/Gq+HQqEQmIiofSyoiNzQrNkTEGK0T1BhlSnw0YErghMRERF1js1mwweni6X2EGMpJk4eIzARUcdYUBG5IZVSiUcGK6X2AU0Mrp3m0D8iInJ+Zw4cx1mdo3dqyehAyOX8yErOi3cnkZuaPm0cokyOxX0/OFEAm80mMBEREVHHrGYzVl1ukNoJphKMSxomMBHRnbGgInJTSoUc3xnhLbVP6qJw7tAJgYmIiIg6dnjXYVzThUjt76RGQSaTCUxEdGcsqIjc2MSJIzHMVCq1V1+qhdVsFpiIiIiobeYGAz7Mc0w8kWItwvDhg8QFIuokFlREbkwmk2FpcrjUztaFIWPvUYGJiIiI2rZ72yHkawMAAHKbFUumxwtORNQ5LKiI3FxCQizGWxyzJX2Qa4O5sVFgIiIiopYM1TX4pNJLaqcpShEdEyYwEVHnsaAiGgC+M3UIZDYrAOCWNgj7th8SnIiIiMhh85bDKFfbn/tVWc341qwEwYmIOo8FFdEAMDg2AtPkjmepPi7Vw1hTKzARERGRXU1RMdYagqT23boKBAf7C0xE1DUsqIgGiEdnjoTSap+QolTji61bDwtOREREBKzbkYlalR4AoLcY8NDcJMGJiLqGBRXRABEWFoi5Gse6VJ/XB6CupLSDM4iIiPpWWU4uNtoipPZ9vg3w8fYQmIio61hQEQ0gD88dC63FCACoVnlg/fbjghMREdFA9uneczAq1AAAH3Md7pkzXnAioq5jQUU0gPj5euEeb8ezU+st4ai4dVtgIiIiGqjyz1/CLmWM1H44zAa9Ti0wEVH3sKAiGmDum5sEL3MDAMCg1OCTXVmCExER0UBjs9mw6shNWOT2hXxDTNWYO3Oc4FRE3cOCimiA8dRr8VCISWrvUEXj9sUrAhMREdFAk51xBkd0jt6pR+O0UKsUAhMRdR8LKqIBaMGsJISYqgAAVpkCqw7niA1EREQDhtViwX+zHJMkDTGVYfqURIGJiHqGBRXRAKRWKfGdIRqpfUwbhfPHzwpMREREA8XRfcdxSRcmtR8fFwS5TCYwEVHPsKAiGqCmTBmNOKNj2vT3zlbAarUKTERERO7OZDRi1U2b1E42F2L0mHiBiYh6jgUV0QAll8vx3bGBUvuyNgSH958QmIiIiNzd9h3HUKDxAwDIbVY8NiVWcCKinlOKDtBcdnY21q5di8uXL8NsNiMyMhLz5s1DWlpat66XmZmJDRs2ICcnBwAwaNAgLF68GOPHt73GgdVqxc6dO7F//37cvn0bJpMJfn5+SEhIwP3334+wsLA2zyNyVQnjhmPCmT04rgoHAKzOsWJio+kOZxEREXVdXU0dPi3VAyp7exYKEB03S2wool7gND1UGRkZWL58OU6fPo2YmBiMHTsWhYWFWLlyJd5///0uX2/Lli14/fXXcfnyZQwbNgyjRo3CtWvX8Prrr2PLli2tjrfZbPjTn/6Ed999F7du3cLw4cMxYcIEKBQK7Nu3D88//zyuXbvWG2+VyKk8NnUI5DYLAKBQ7YttO44JTkRERO5o7bbjqFZ5AAC0lkZ8azYnoiD34BQ9VLW1tVi5ciWsViuee+45pKSkAAAqKyvx0ksvYfPmzRg/fjwSEhI6db38/HysXr0aKpUKy5cvR3x8vLT/xRdfxOrVqzFu3LgWPU6ZmZk4ceIEgoOD8bvf/Q6+vr4A7L1Wq1evxubNm7Fq1Sq8/PLLvfvmiQSLiovBnMM7sV0RBQBYU6bHQ2UV8AnwE5yMiIjcRWFuHjY0+ANfzYx+r6YUAaFjxIYi6iVO0UO1Z88e1NfXIzk5WSqmAMDX1xdLliwBAGzatKnT19uyZQssFgvmzJkjFVMAEB4ejvvvvx8WiwVbt25tcc6FCxcAALNnz5aKKcD+nMmDDz4IAOyhIrf1yNyx0FoaAQA1Kj3+++luwYmIiMid/O2LgzAq1AAAX1Mt7luQcocziFyHUxRUmZmZAIDU1NRWryUlJUGlUiErKwtGo7FT1zt58mS715s0aVKLn9lEpVK1ez3ZV1N5enp6durnE7ka/9Ag3OfhWBNkncEfeVeuC0xERETu4uKJM9iFUKn9SLARei9+piL34RQFVW5uLgAgNrb1TC9KpRLR0dEwmUzIz8+/47Xq6upQWmqfCnrQoEGtXg8ICICXlxdKSkpQX18v7R89ejQAYPfu3aisrJT2W61WfPbZZwCAGTNmdPo9Ebma++5OgZ+pFgBgkqvw9/VHBSciIiJXZ7PZ8NddF2GT2T9yRhrLMWfuRMGpiHqX8Geo6uvrUVdXBwDw9/dv8xh/f39cu3YNpaWlbRZJzTUVUx4eHtBqtW0eExAQgJqaGpSWliI6OhoAMGrUKCxatAibNm3CT37yE4wYMQJarRY3btxAeXk5FixYgIcffrjT7+vZZ59ttU+tVuO1114DAAQGBrZ6vT8plfa/+qCgIKE5yLl8N+4qVty0b+9RROAb2dcxeiqHZZADf3dQe3hvUFvSN+3GSW2E1H4yOZSzJlML7vC7Q3gPlcFgkLY1Gk2bxzTtb37sna7X3rU6ut7SpUuxdOlSmM1mnD59GkePHkVRURHCw8MxcuRIKBSKO/58Ild23+I0xJjsQ/9sMjneOXgTVrNZcCoiInJFpoYG/PVctdROtJYhbSZ7p8j99EoP1Ztvvolbt2516ZynnnoKcXFxvfHjW7DZ7KtvNz331Fkmkwl/+ctfcOzYMTzwwANIS0uDl5cXrl27hv/+97/405/+hO9973uYP39+p663YsWKDl8vLS2VsorQ9C1ASUmJsAzknB4fE4CXL1gBABe0Ydjw8XpMmT9dcCpyFvzdQe3hvUFft/mLXbiljQQAyGxWfHdSjDSSiKiJ6N8dMpmsx72mvVJQlZSUdOr5puYaG+0zijUfltfY2Ai9Xt+pY9uj0+kAdNyb1db11q1bhyNHjrQa2peQkIBf/epX+OlPf4qPPvoIU6dO5eQU5NaSxsVjwoWDOA77sNT381VIrqmFhg8QExFRJ9WUlODjGj9pEd+5mgoMiRspNhRRH+mVgqrpuaDu0Ov10Ov1qK+vR3l5eZsFVXl5OYDOPXfUdExdXR0MBkObRVhZWVmr6x04cABA2zMDBgYGIj4+HllZWbh27RrGjOG6CeTenlmcjMe+vAaLXIEijR82bT6EBx+ZJzoWERG5iDVbTqBGNRgAoLUY8eQ3ZgDgEHJyT8KfoQKAmJgYAMD1662naTabzcjNzYVKpUJ4ePgdr+Xh4SEVSjk5Oa1eLysrQ01NDQIDA1sUb01FVlsFHeDo+aqtrb1jBiJXFxsbhXs8HePeP2sMQcWt2wITERGRq8i7eBlb5FFS+9FQM4KDuFg8uS+nKKiSkpIAAEePtp6m+eTJkzCZTEhISIBare7x9Y4cOdLimCZNi/m2tXiv1WrFjRs3ALj2DCREXfGjb94FT3MDAKBBqcWHu7IEJyIiImdns9nw3qEbMMvtg6ACTTVY8tBMwamI+pZTFFSzZs2CTqfDiRMncOzYMWl/VVUVPvjgAwDAokWLWp23bNkyLFu2TBoS2GTBggWQy+XYuXMnLl++LO0vKCjAunXrIJfLsWDBghbnTJgwAQCwZs2aFs+DWa1WfPTRRygpKUFQUBCGDBnS8zdM5AJ8fDzxzVCT1N6tisH1kyyqiIiofWfSM5Chi5Hajw/VQqvp3BfiRK5K+DpUAODp6Yknn3wSb731FlasWIGRI0fCy8sLWVlZqKurw913343ExMRW5zUVPuavTescHh6OJUuWYNWqVVi+fDlGjx4NhUKBs2fPwmg0YunSpa2GDz700EM4c+YM8vPz8fOf/xzx8fHw9PRETk4OioqKoFar8eSTT3LqdBpQ7p6djG2rM5Cn8oVVJsd/MvPxypgRkCuc4lcHERE5EXOjAf+50gjYn5LAMHMppk6eIjYUUT9wmk9FqampePnll7F27VpcuXIFZrMZERERmDdvHmbO7HpX8aJFixAaGoqNGzfi4sWLAIDY2FgsXrwYycnJrY738vLCq6++io0bN+L48eO4evUqzGYz/Pz8MGPGDNx7772IjIzs8fskciUqhRzfTfTF7y7Z21n6SGTsOoLUedPEBiMiIqeza/NB3NQ5Pit9f1JMl5exIXJFMpvIxZAGsIKCAq5DRU6r+f1hs9nw/1YfxGmFfV+YoRzvPDIGag8PkRFJEP7uoPbw3hjY6kpK8eSmHFSp7UtszJAX49lv2dcw5L1BHRF9f/TGOlRO8QwVETkvmUyG780YCrnNvthvgdYfmzcdFJyKiIicyWdbM6RiSm0x4TtzRwtORNR/WFAR0R3FxIRirrpMaq8xhKDidtcW8yYiIveUd+kKNsqipfYDvrUICvAVF4ion7GgIqJOeXT+eHiYDQCAeqUWq3dyxj8iooHOZrPh3UM50jTpAaYa3Dev9bPqRO6MBRURdYqPtx7fCjVK7d3qGGRnnhOYiIiIRDuRfhyZWscivt8dqoFOoxKYiKj/saAiok67e9Z4RBsd677963QZLBZzB2cQEZG7Mhoa8e41i9QeZSrB1Mmtl7khcncsqIio05RKBX4wxl9qX9GGYM/ODIGJiIhIlPVbjqBA4wcAkNus+MEUTpNOAxMLKiLqktFJwzHF7JiQYnWhBrU1tQITERFRfystLMXntX5Se768EIOHRndwBpH7YkFFRF32+OyRUFvsz1NVqTzwyebjghMREVF/en/nWRgUGgCAt6kOj97NiSho4GJBRURdFhwRiod0pVJ7szkYuVdyxAUiIqJ+c+7kBaTLw6X2twPr4eXnLTARkVgsqIioW+5bOAkhxkoAgFWmwL8O3IDVahUbioiI+pTZaMK/TjsmJ4o1lGD2vEkCExGJx4KKiLpFo9fhe0McU+Oe1YThyJ5jAhMREVFf27H1EHI0gVL7h+ODoFQpBSYiEo8FFRF1W8q0JIw1FUrt/+TKYaiqFpiIiIj6SlVhMT6scAztS7PmY8TY4QITETkHFlRE1G0ymQw/mBkPhdW+DkmpxgefbzosOBUREfWFD7dmolalBwDoLI1Yevc4wYmInAMLKiLqkciYcCzSOcbTr7NGIu/CZYGJiIiot2VnnMIOVYzUfjigAQGBfh2cQTRwsKAioh57ZOFE+Jvsa1GZ5Ur889BNWC1mwamIiKg3mBsN+PvZathk9o+NUaYKLJo7QXAqIufBgoqIekyv0+D7wzRS+7Q+Coe2HxSYiIiIesu2Dem4rguR2j9KDoZapRCYiMi5sKAiol4xJTUBYy0lUvs/RXrUl5UJTERERD1VnnsbH9YHSe00eTESE4YITETkfFhQEVGvkMlk+NGc4VBZ7UP9ytXe+HgTp1EnInJVNpsN7+3MQr1SBwDwMBvwOCeiIGqFBRUR9ZrwsCA84F0jtTcpYnDjZJbARERE1F1Z6RnYrx0stZdEmOHn6yUwEZFzYkFFRL3qgbuTEWqqAgBYZQr8LbMEFpNRcCoiIuoKY10d/nHFJLWHmMowd2aSwEREzosFFRH1Kq1ahR+McUylm60Px+5N6QITERFRV63fcAC3dYEAAJnNiienREOp4MdGorbwXwYR9brkcfGYBMcEFauq/FCVXyAwERERdVbR5WtYYwqX2vPUZRg6NEpgIiLnxoKKiPrE9+YlQmuxD/WrUXlg9bZTghMREdGd2KxW/Hv/VRgVagCAj7kOSxYmC05F5NxYUBFRnwgO9MUjwQapvVMTi4uHTwhMREREd5Kx8xAy9DFS+7uxKnh56AQmInJ+LKiIqM8smpOMaFOF1P7rxQYYGxoEJiIiovbUV1TiH/lqqT3KUoq0qYkCExG5BhZURNRnVAo5fpwaJrVvaYOwbsNBgYmIiKg9H2w8hjK1DwBAaTXjf2YOhUwmE5yKyPmxoCKiPjVi+CDMVxRJ7c+Mobh9NUdcICIiaiX79AVskTsmnnjIqxLRUSECExG5DhZURNTnvrNoIvxMtQAAk1yFv6XnwGq1Ck5FREQAYDKa8NeT5bDJ7B8LI4wVeHBBiuBURK6DBRUR9TlPTx1+GOsYNnJOE4rduzIEJiIioibrtxzBTU2g1P7xWD+o1SqBiYhcCwsqIuoXk6YlYaIxT2q/V6BBRWm5wERERJR/Mw+f1jgWY59rvY2EccMFJiJyPSyoiKhfyGQy/HDuKGgtjQCAWqUO/9nMtamIiESxWq342+7LMMrtvVG+xhosXcQ1p4i6igUVEfWboIhQfMe3UmqnKyOQeYBrUxERibBvx2Gc1ThmYn0i2gIvP19xgYhcFAsqIupX8xdMwdDGYqn99ytmNFRVC0xERDTwVBUW4z+FeqmdbMzHlJkTBSYicl0sqIioXymVSvzvtBgorBYAQLHGFx+vPyQ4FRHRwGGz2fDulpOoUdkLKq3FiB/NT4Bczo+FRN3BfzlE1O8GD43BfR6OCSk2ymNw+RifpyIi6g8ndx/Gfs0gqf1oUD2Cw4LFBSJycSyoiEiIhxemIsxUBQCwyuT4c1YtjLV1glMREbm3upJSrMxVSu0hpjIsnMuhfkQ9wYKKiITQalT4cXKQ1M7VBeGzL/cLTERE5N5sNhve35iBUo0PAEBpNeMnabFQKvhxkKgn+C+IiIQZnRCL+eoSqf2FLQbXM88KTERE5L7OpB/Ddk2s1H7ItwaDB4V1cAYRdQYLKiIS6rFFKQgy2Wf5s8gVeOdUBUz19YJTERG5l4aKCqy8ZpPaMaYKPHh3isBERO6DBRURCaXXqfG/4/yk9g1dCNZx6B8RUa/6YP0RFGnsv2vlNguenhYFtZIfA4l6A/8lEZFw48YMxWylY+jfp5Yo5J45LzAREZH7uHAwA5uVg6T2/Z5ViBsSKS4QkZthQUVETuG7iybA31QLADDLlXgnoxhmg0FwKiIi12aorsafs02wyewf+SJNlfjmQs7qR9SbWFARkVPw9NDix4leUvuKPgwbv9wnLhARkRv4eN1B5GsDAABymxVPTwqHRqW8w1lE1BUsqIjIaUwYPwxpcsfQv4+MEcg7d0lgIiIi15V9NBMbFIOk9j36CgwbFi0uEJGbYkFFRE7l+/eMh6/ZvsCvUaHCn4/kw9xoFJyKiMi1NNbU4s/nGmD9aqhfmKkKjy7iUD+ivsCCioicirenHv8zQie1L+rDsWE9Z/0jIuqKD788iFu6QKn91IRgaNUqgYmI3BcLKiJyOpMmjsQ0WbHU/rAxDDcvXBWYiIjIdZzLOIMNMsfQvoXqUiSMGiwwEZF7Y0FFRE7ph/eMh1+zWf/eOlIAI4f+ERF1qL6mFu+cb5Bm9YswVmDpPRzqR9SXWFARkVPy9vLAUyM0UvuGNghrNhwWmIiIyPn9d0MGitS+AOwL+D4zMRharVpsKCI3x4KKiJxWckoi5trypPYXjUG4fO6KwERERM7rxNEs7EC41H5QXYxho4YITEQ0MLCgIiKn9t17UxDSWAEAsMoUeDujBIb6BsGpiIicS3V5Jf5yybEY+mBDCR6+d4rAREQDBwsqInJqei9PPJ3oCZnNCgDI0/hj9dqDglMRETkPm82Gf6w/jgqVfXF0pdWMZZPDoNZwqB9Rf2BBRUROL2FCIhYrC6X2JlkUzh46ITAREZHzOLj9EA6qo6T2o15lGDQiTmAiooGFBRURuYRv3z8FUcZyqf1/2WbUlZQITEREJF7ZzVz8vVAvtUcYi3HvoqkCExENPCyoiMglaDQaLJsSAYXVAgAo1fjiX+uPw2a1Ck5GRCSG1WTEX3ZcQq3KXlBpLUY8PXc4lEqF4GREAwsLKiJyGXHxMXjYt0Zq79XF4sCmvQITERGJs2ntHpzUOxbw/W6kGeERwQITEQ1MLKiIyKU8tGAihpnLpPbfKvxRfPmawERERP3vxqksrDJGSO1kWwnmzRwnMBHRwMWCiohcilIhx0/nj4DO0ggAqFfq8Nb+HJgNhjucSUTkHgzV1ViRWQmTXAUA8DXX4al7kiCTyQQnIxqYWFARkcsJC/HHD2MdHxwu6CPwxRd7BCYiIuofNpsN739+ALm6IGnf02O84efjITAV0cDGgoqIXNLMKYmYJi+V2p/YYpB9JFNgIiKivndi10Fs0QyR2ou0ZRg/dqjARETEgoqIXJJMJsP/LE5GkKkaAGCVKbDiQiPqy8rucCYRkWuqyL2NP9/WSe0YUwWWLk4RmIiIABZUROTCPD20+GlKMOQ2+9TphVp//PPLY5xKnYjcjtVkwjvbL6BK7QkAUFtNeO6uWGhUSsHJiIgFFRG5tFEjBuFB7yqpvVcbi4Nb9okLRETUBzavazlF+mNhRsREhwhMRERNnOprjezsbKxduxaXL1+G2WxGZGQk5s2bh7S0tC5dp7q6GsePH8fVq1dx9epV3Lp1C1arFc888wymTJnS4bm3b9/GmjVrcP78eRgMBoSGhmLmzJlYsGAB5HLWn0TO6JGFKTjzwRFcVgYAAFaW+SH+6nWExMUKTkZE1HM3Tp3D+43h0tfg420lWDhrqthQRCRxmgohIyMDy5cvx+nTpxETE4OxY8eisLAQK1euxPvvv9+la126dAn/+Mc/sHv3bty8eRPWTg7/uXz5Mn71q1/h6NGjCAkJQXJyMmpqarBq1Sq89dZbsNls3XlrRNTHlAo5np03vMVU6n/adxMmTqVORC6uoaoab5ysajFF+k8WjeMU6UROxCl6qGpra7Fy5UpYrVY899xzSEmxP2BZWVmJl156CZs3b8b48eORkJDQqev5+vpi7ty5iIuLw5AhQ7B+/Xqkp6d3eI7FYsGf//xnNDY2YunSpVi0aBEAwGAw4He/+x2OHTuGffv2YebMmT17s0TUJ8JCA/CjwXl4O9fevqwLw4dfpOPxb88VG4yIqJtsNhv+se4o8rSOoX5Pj/aCn6+nwFRE9HVO0UO1Z88e1NfXIzk5WSqmAHthtGTJEgDApk2bOn29+Ph4PPHEE0hLS0NUVFSnvsXJyMhAUVERYmJipGIKALRaLb7//e93OQMR9b+Z00bjLlmR1F6HaBw/eFJgIiKi7tu94yj2qhzF1P3aEowfFy8wERG1xSkKqsxM+9oxqamprV5LSkqCSqVCVlYWjEajkAyDBw9GSEgIbt26heLi4j7LQEQ998P7UxHVWC61/+8aUFJY2sEZRETOJ/f6LfyzyLFY7zBjMb597ySBiYioPU5RUOXm2sfoxMa2foBcqVQiOjoaJpMJ+fn5fZbh5s2bAOzFU1ua9jcdR0TOSafT4OfTwqG22L+AqVHq8ebW8zCbzIKTERF1jqGuAW/sy0WjQg0A8DQ34Gdz46FSOsWTGkT0NcL/ZdbX16Ourg4A4O/v3+Yx/v7+uHbtGkpLSzFo0KA+yVFaav8GOyAgoN0MzY+7k2effbbVPrVajddeew0AEBgY2J2YvUb51S/loKAgoTnIObn6/REUFISf5G7Bm4X2DyMX1SH4bP0hPP3kQ4KTuT5Xvzeo7/De6B02mw2vfLgGuZpwad8vR6gwamyiwFQ9w3uDOuIO94fwHipDs1m4NBpNm8c07Tf04YxdTdduL4NWq+3zDETUe+5/eD5mWh292p8agnBg6z5xgYiIOmHTZ1uxHY5i6gFVIe5ayAmxiJxZr/RQvfnmm7h161aXznnqqacQFxfXGz++X3R1yvQVK1Z0+HppaanQadibvgUoKSkRloGcl7vcHz+8PxXZn51BvtoPNpkcvz9Xi7eDMhEQE33nk6lN7nJvUO/jvdFzty9cxlu5aunT2VBjCb710BSX/zPlvUEdEX1/yGQyhIWF9egavVJQlZSUdPn5psZG+3oxTT0/Tfv0en2nju1tWq0WdXV10s8SkYGIepfeU4+fT43AL45WwyRXolrliRU7LuPlJUFQ6nSi4xERSQyVVXjjcAEMuhAAgIe5AT+bNxxqjUpwMiK6k14pqJqeC+oOvV4PvV6P+vp6lJeXt1lQlZfbZ+zqy+eOAgMDUVdXh7KyMsTExAjJQES9L3ZoNJ7IzcTfCu2/7s7pI/Hhp7ux9LGFXBiTiJyC1WLGPz4/iBzdEGnfT4arERruus+UEA0kwp+hAiAVMNevX2/1mtlsRm5uLlQqFcLDw1u93tsZbty40ebrTfvbKraIyLnNuysJU+WOCWXWquJwZMs+cYGIiJrZsXYX9jQrphbqyjEpZZTARETUFU5RUCUlJQEAjh492uq1kydPwmQyISEhAWq1WkiGGzduoKioCJGRkQgODu6zDETUN2QyGf73/hREmiqlfe+U+uF21gVxoYiIAFw+cgL/MkRK7WGWMjy+uPWamETkvJyioJo1axZ0Oh1OnDiBY8eOSfurqqrwwQcfAAAWLVrU6rxly5Zh2bJl0nC8npg4cSKCg4Nx8+ZNbNq0SdpvMBjw7rvvtpuBiFyDXqvCL+cMgc5ifx6yQanFa8fK0NALvz+IiLqj6nYeXr9ogVluH5Lsa67D8/ckQq10io9nRNRJwtehAgBPT088+eSTeOutt7BixQqMHDkSXl5eyMrKQl1dHe6++24kJrZef6FpIgyzufWCnb/+9a+l7cLCQgDAp59+ii1btgCwL9T7xBNPSMcolUr85Cc/wW9/+1usWrUKR44cQWBgIC5duoSKigpMmDABaWlpvfm2iaifRUUE4enhpXj9ir19SxeEv649hmcfnwM5F8wkon5kNhjw5tYLKNVHAQDkNgt+luyPAD9vwcmIqKuc5hNEamoqXn75ZaxduxZXrlyB2WxGREQE5s2bh5kzu77+wpUrV1rtKywslIorlar1rDnDhg3Dq6++ijVr1uDChQvIyclBSEgIFi1ahIULF0Iu5zdGRK5u8sQRuD//CNbV+QEADugGI/7znVj8yN2CkxHRQGGz2fDxp7txRu94bmppUAMSR/G5KSJXJLOJXAxpACsoKOA6VOS03P3+MFusWP7hYZxT2GftVFgt+O3gOoyaOlFwMufn7vcGdR/vjc47um0/Xi0LkdqTUILnH53qtjOP8t6gjoi+P3pjHSp2uRDRgKNUyPGze8fB31QLALDIFXjjClB+M1dwMiJyd/kXLuP/inykdoSpEk8/kOK2xRTRQMCCiogGJD8fDzw/KRhKqwUAUKH2xh93XIWxoUFwMiJyV/UVlXjtcCHqlVoAgNbSiF/OioVe13ezGBNR32NBRUQD1vBh0fh+uEFqX9SH41+fpsNqtQpMRUTuyGIx4+21x3FT51h+5el4JaKjuBwLkatjQUVEA9rddyXhLnmx1N6hisHWzQcFJiIid/TpF+k4po2S2vfpyjAllZNQELkDFlRENKDJZDL8z4OTEG90PAz7bmUAzmZy0V8i6h2H0k/iU1O41E4yF+I7904SmIiIehMLKiIa8DRqFX65cAT8jdUA7JNU/PGcAYU38wQnIyJXd/3SNfxfjmOVmojGcjz3QDKUCn4EI3IX/NdMRAQgIDgQvxrnBZXVBACoUerxh903UF9dLTgZEbmqyqJS/OFIKRoV9kkn9OYGvDA9HJ5enoKTEVFvYkFFRPSV+LEj8L8htVL7piYQb392FBaTSWAqInJFxoYGvL7xLErU9inS5TYrfh5nQ2TcILHBiKjXsaAiImpm5txJuF9VILWPaaPxySc7hS7ETUSuxWa14N8f78UFneO5qaXe5UiakiQwFRH1FRZURERfs+TBGRhvLpLaa+SxOLhhl8BERORKtq7Zhu2aWKmdhiLcd88UgYmIqC+xoCIi+hqlQo5nH0pBhKlS2vdOVSguHTwmLhQRuYQT2/fjX+ZBUnuoqQw/fngyZDKZuFBE1KdYUBERtcHTQ4sX5g2Fp9m+8K9RocLvryhQeDFbcDIiclbXj5/CG4W+sMoUAAB/cy1+dW8iNCqV4GRE1JdYUBERtSMyLADPT/CF0moBAFSrPfHK4RLUFBYKTkZEzqbseg5+l2WEQakBAGgtRvxmRiQC/LwFJyOivsaCioioA6NHDsL/DrZK7TxtIF7fmAVjbW0HZxHRQFJfVobf7rmJMo1jRr+fJeowZFCo4GRE1B9YUBER3cFdUxPxsI9jPaosfRRWfpIOK6dTJxrwzAYD3lx7Ajd0IdK+J8IaMWHcUIGpiKg/saAiIuqERxdOwHRFqdTeq4vFZ5/s4HTqRAOYzWrFfz7ajRP6GGnfPfpyLJw1TmAqIupvLKiIiDpBJpPhJw+mYqTZUVR9JB+C9I17BKYiIpE2rtmBzZohUnsiSvH44lSBiYhIBBZURESdpFYp8csHkhBmqpL2vVMVgnMHjwtMRUQiHNmejv+ao6X2EFM5nn0oBUoFP1oRDTT8V09E1AU+Xnq8NG8IvMz1AACzXInfX1PixvnLgpMRUX85l3EGK4r8YJXZP0YFmmrwm3sTodNwenSigYgFFRFRF4WHBeKFCX5QW+2TUtQrdXj5eBUKb3E6dSJ3dyM7B3+4aIVRYS+e9GYDXkyLhL+fl+BkRCQKCyoiom4YOXIwnhtig9xmn1K9QuWFl3fdQGVZheBkRNRXivIK8cqREtQpdQAAldWEX4/VY9CgMMHJiEgkFlRERN2UOnk0/ifAUUDlq/3w2/VZqK/hGlVE7qaqpAwvb7+OcpW9J0pus+K5aCMSxsQLTkZEorGgIiLqgXl3T8GjqjypfVUTjNfXHIOxwSAwFRH1poaqKvx2wznkafylfT/yKcGkGeMFpiIiZ8GCioioh77x0EwskOVL7dPaCPz5o32wGI0CUxFRbzDV1eH1zzJwRetYuPdRbSHm3zNDYCoiciYsqIiIekgul+P7D8/AFKtjUop07SD8d/V2WE0mgcmIqCcsBgP+/NE+nNJFSfvuVhTiG/dPF5iKiJwNCyoiol6gVCqw7JEpGG0pkfZt1A7Fpx9shc1qEZiMiLrDajTinx/sxH69Y+HeySjGE9+YDrmcH5+IyIG/EYiIeolapcIvv5GCWLNjoopP1PFY++EW2Gw2gcmIqCusZjPe+2A7tumGSvsSrGX46cOTuXAvEbXC3wpERL3IQ6fGSw+MRaS5Stq3Sj4Umz9hUUXkCmxWKz75cCvWaxzFVLylHL/+xkSoVUqByYjIWbGgIiLqZX5eOryyeCRCTdXSvn9Zh2DX2u0CUxHRndhsNqz9aAs+VTqKqcHmCrz0UBL0WpXAZETkzFhQERH1gQA/L7yyKB6B5hpp318borF/w26BqYioPTabDZs/3YZVsjhpX6SpEv/vgTHw0msFJiMiZ8eCioioj4QE+uKVubHwNdcBAGwyOd6uDsWRrfvEBiOiVnau3Yl/WQZL7VBTNV65dxR8vfQCUxGRK2BBRUTUhyLCAvDKXZHwMjcAAKwyBf5UGogTuw8LTkZETfZv3IOVDZFSO9BUg1cWxiPAz0tgKiJyFSyoiIj6WExUCF6eFgy9xQAAMMuVeC3fC5npxwUnI6L9Ww/g7aoQ2GT2j0R+plr8dt5ghAT5ig1GRC6DBRURUT8YEhuB5Sl+0FoaAQAmuQp/uKnD8cOnxQYjGsD27DyKt8sCYJUpAADe5nq8clcUwsMCBScjIlfCgoqIqJ8MHxaDl5I8paLKLFfitWtKHDt4UnAyooFn146jeKfIG9aveqa8zfV4ZVoIoqNDBCcjIlfDgoqIqB+NShiC/zdGA53ZMfzv9RsaHN5zVHAyooFj+6Z0/KXYWxrm52OqxW8nB2BwbITgZETkilhQERH1sxFjhuPlRBX0X01UYZEr8Ea+Fw5sTRecjMi92Ww2bF27CyurgqViytdUi99N8sOgoTGC0xGRq2JBRUQkwLCkUXhlvCc8ms3+t6IsAPu/3Ck4GZF7sq8ztRV/bzabn7+pBr+fFozoYUMEJiMiV6cUHYC6xmaz9ep1eut65F766v6QyWS9ej1XNzRhKF5R5WD50QrUKnWwyhR4uzYC5jWbcdc3FvDPi6iX2KwWrP9gE/6rGCbtCzDV4Hd3RSE8OlRgMiJyByyoXIDZbEZDQwOMRmOvfcCtqKgAAFit1l65HrmXvro/ZDIZ1Go1dDodlEr++gGAuGGD8FuVCi8dKEaNUgerTI53TENQ++FGLH50IWRyheiIRC7NajLiw9Vb8bnGUUwFmmrwu3mDEcbZ/IioF/ATjZMzm82oqqqCVquFr68v5PLeGaXZ9GHWbDb3yvXIvfTV/WG1WmEwGFBVVQUfHx8WVV+JjY3A71RqvLTnFqqUegDAf2TxqHl/Mx5dcjfkKpXghESuyWxowD8/2IXtOkcxFWyuwW8XDkUo15kiol7CTzNOrqGhAVqtFh4eHr163aahRBxSRG3pq/tDoVBI93JDQwO8vLx69fqubFBUEF6dr8TybVdRorT/uXymjkf1qu344bdnQ6nVCk5I5FqMtTV4++MDOKQfKu2LMlfi5XsTEeDbu/9PJaKBjZNSODmj0QgtP0iRm9FoNDAajaJjOJ2IED+8du8IRJkrpX3btXFY8cFeGOvqxAUjcjEN5eX4/cdHcEgfK+2Lt5TjDw+NZTFFRL2OBZUTs9lssNlsvTbMj8hZKBQK6f6mlgJ9PfH7B8diqLlc2ndINxi//+gwGr56to2I2lddUISXvjiD0/poad9YWxleeWQCvD34BSUR9T5+UicicjI+nlq88s0JGGMtlfad1kdh+eenUVVYLDAZkXMrzcnFrzdfxmV9mLRvirwUv35kEnRqPotIRH2DBRURkRPSa1X4zSOpmCxzFFXZ+jD8cstV5N+4LTAZkXO6fuEqfr63ELm6IGnfXHUZnn14MtRKftwhor7D3zBERE5KrVLiuW9Oxly1o6jK1/jj+f3FuHj2ssBkRM4l8+hZ/OpEHcrV3tK+h7wq8OOHJkOp4EcdIupb/C1DROTElAo5fvzQFHzLs0zaV63S48UzjTiw+6jAZETi2Ww2bNuUjt9dVcCg0AAA5DYrfhhYhe8snsSZbImoX7CgIpdVX1+Pf/7zn3jooYcwZswYDBo0CCNHjsQ999yDN954A3l5eaIjDhgpKSmIiIgQHcNtyWQyPHLvFDwTVAGl1b42mEmuwp8KffHFJ9thtXA9ORp4LCYj3v9gO/5WFQyrzL4AttbSiBcGGbBwXorgdEQ0kLCgIpeUmZmJqVOn4uWXX8bp06cxbNgwLFy4EOPHj8fNmzfx9ttvY9q0aUhPTxcdtd8dPnwYERERWLZsmego1MvumjsJ/2+4DZ7mBmnfKksM/vbedphrawUmI+pfhooK/Om93VgnHyTt8zPV4A9JOkyYmiQuGBENSFzYl1zOhQsX8PDDD8NgMOB///d/sWzZMuj1eul1q9WKbdu24fe//z0KCgoEJiXqfYkTEvFaUB5+u/c2itQ+AIAd2iEo/vAQfrZgFLyiIgUnJOpbFVev4tU9N5HtMVjaF2Msx2/ujkdwaKDAZEQ0ULGHilyKzWbD008/DYPBgOeeew4vvPBCi2IKAORyORYsWICtW7dizJgxgpIS9Z2oQRH44wMjEW92PFd12jMGv9ieg9zM0+KCEfWx7ANH8Wx6GbI9HEOMx1lL8eo3x7OYIiJhWFCRS9m3bx8uXryIsLAwPP300x0e6+3tjeHDh0vthoYGvPXWW7jrrrswZMgQDB8+HA888ADWr1/f5vnNnwt67733pPNSU1OxcuVKaVHarKwsLF26FKNGjUJ8fDy+973v4fbt1tNaL1u2DBERETh8+DD27NmD++67D0OHDsXIkSPxxBNP4OrVq63OefPNNxEREYFPP/30jhmbfsY3vvENAMBnn32GiIgI6b8333yzxbm3bt3CL37xC6SkpGDw4MFITEzED37wA1y4cKHNn2U2m/HnP/8ZU6ZMQWxsLCZNmoQ//vGPMBqNbR5PfcvXxwu//VYKJikci/3m6wLxi3PA0U27uGgyuRWb1Ypdn23Fr3M8UK7xkfbP1VXi149OhodeIzAdEQ10HPJHLmX37t0AgEWLFkGp7PztW1tbi2984xs4e/YsAgICMGvWLDQ0NODQoUM4duwYMjMz8corr7R57vLly/HBBx9g3LhxiIqKwtGjR/H73/8e9fX1mDFjBr71rW8hKioKkydPxoULF7B9+3ZkZ2dj165d0Ol0ra63adMmrFq1CmPGjMGcOXNw8eJFbN26FYcOHcLnn3+OUaNGde8PB8DEiRNRUlKCffv2YdCgQZgwYYL0WvPrZmRkYOnSpaipqcGwYcMwZ84cFBYWYuvWrdizZw8+/PBDTJ06tcW1f/zjH2Pz5s3w8PBAWloabDYb/vnPf+LcuXP88C6IVq3EL76ZijXbT+LjMg8AQINSi1erIvHIqs14+FtzoVCrBack6hlTfR3++8lebNbESV8Dy21WfD/KjIXTUziTHxEJx4LKRdlsNqChrvvnf1WM2Mz9PDuYzqNH//M7d+4cACAxMbFL57322ms4e/Yspk2bhnfffRceHvYPn1evXsWDDz6Id999FzNmzMCsWbNanbtp0yZs2bIFw4YNk86ZO3cu/v73v+Pzzz/H888/jx/84AcAAKPRiCVLluDQoUPYsGEDvvnNb7a63vvvv48//vGP+Pa3vw3A/nf56quv4q9//Suee+45bNu2rUvvrblHH30UgwYNwr59+zBhwgS8/fbbrY6pqanBj370IxgMBvzjH//AokWLpNfS09Px2GOP4amnnkJGRgbkcvunly+//BKbN29GTEwMvvjiC4SFhQEAcnNz8cADD/BZNYHkMhkemT8eg09dxVtZdWj4auroT5RxuP7+Xiy7bzw8gjgUilxTZV4+3thyAef0cdI+b3MDfj4xAKNHRAtMRkTkwILKVTXUwfrMo90+XdQgLfn/fQToPbt9fkWFfXhTQEBAp8+pr6/Hxx9/DLlcjj/84Q9SMQUAcXFxeOaZZ/Diiy/iP//5T5sF1S9+8QupmGo6Z9asWdiyZQsiIiKkYgoA1Go1nnjiCRw6dAhHjhxps6BKTk6WiinAPiX2z3/+c6xbtw5ZWVk4ceIEkpOTO/3+uuqTTz5BcXExnnrqqRbFFABMnz4dS5cuxb///W/s3LkT8+bNAwCsWrUKAPDzn/9cKqYAIDo6GsuWLcPzzz/fZ3mpc1LGxeGPwaX4w65rKFDah0Rl6GPwiw2X8KtJwYgcGS84IVHXXMs8i1fP1KNE75hoZbC5Ar9aNAohAd4dnElE1L/4DBW5lO4MLTt79iwMBgPGjh2L2NjYVq8/+OCDAIDjx4+3ef1p06a12hcdHd3uazExMQCA4uLiNvPce++9rfapVCosWLBAytGXmqaSnz9/fpuvT5w4EQBw6tQpAIDJZMKpU6cgl8uxcOHCVsffd999fROUuiw6IhBvfGMMxttKpX23tYH42fF6HNjFRYDJNVitVmzdmI7nL8hQovGV9k9XlOK1byWzmCIip8MeKnIp/v7+uHbtGsrKyu588FeKiooAAFFRUW2+7uPjA29vb1RXV6Ompgbe3i3/Z928R6ZJ08yCHb3W3mQNkZFtT2vdlK+wsLDN13tL04QZX++d+rry8nIA9l5Bo9GIkJAQqNt4HsfT0xM+Pj6oqqrq/bDUZV56LV741mR8vOkYPq/1A2B/rupPRVqc/XAPvv/AJGjbeLaPyBnUVddi5ZcZOKgIb/G81NLAOtw3bwqflyIip8SCylXpPOzD57qpaUIHs4BnqHpi1KhROH78OLKysqSepc7qzP+I2zqmv/4H3p3eN6vV2uVzLBYLAHtB1dakGYB96vmkpKQWufhBxnUoFXJ8595JiD2chb9csaBeqQUA7EA4sj/KxC+mhCJyeNwdrkLUv66eOo8/napGgSZc2udtqsOzY70wbuyEDs4kIhKLBZWLkslkPXoWSfZVQSXr74Kqh2bNmoX33nsPmzZtwm9+85tOzfQXEhICwD6BQluqq6tRXV0NvV4PT8/u/5l2VltTqgNAXl4eACA0NFTap1KpANifA/s6i8WCkpKSLv/8sLAwXLt2Dc888wxGjhzZ5jHNC25/f3+o1WoUFxfDaDS26qWqra1l75STmjI5EYMj8vCnXddwTRsMALipDcRzGXX4n0u7kLZ4JmRyheCUNNBZTUZsWbcH/22MhFnjJ+1PMBTg2UWJCAjhpCpE5Nz4DBW5lJkzZ2LYsGEoKCjAO++80+GxNTU1yM7OxujRo6HVanH69Glcv3691XFr164FYH92qD96YTZs2NBqn9lsxpYtWwCgxYQUTcVgW7kPHToEk8nUan9TEdbUE/V1Tc99bd++vVN5VSoVxo4dC6vVKmVsrr11vMg5hMdE4LVvp2CRPF/aZ1Bo8HZdJP78780wFBcJTEcDXU1uLl7/zy78yzQIZvlXX/TZrPimphAvL53GYoqIXAILKnIpMpkM77zzDrRaLd588028+uqrrXpvbDYbduzYgbvvvhunT5+GXq/HI488AqvVil//+tctjr927Rr+7//+DwDw3e9+t1/ew/Hjx/HJJ5+0yPunP/0JeXl5GDlyZIu1o1JTUwHYi75bt25J+2/evInf/OY3bV6/qYfr2rVrbb6+ZMkSBAQE4M9//jM+/fTTVkMN6+vrsWbNGuTn57c4BwDeeOMN6Zk0wN7b1tbU7ORc1FoNfvCtu/DLWBM8zAZp/26PePx0wxVc2nuQa4lRv7LZbDi1dTee2VWAo56OyYJ8TXV4OUGBRx9Kg1LFQTRE5Br424pcTkJCAj755BP84Ac/wF/+8he8++67GD9+PIKCglBdXY2zZ8+ipKQEWq0W4eH2sfi/+tWvcPLkSaSnp2PSpElITU1FfX09Dh8+DIPBgO9///uYPXt2v+RfunQpfvazn+GDDz5ATEwMLl68iOzsbHh6euKtt95qcWxMTAweeughfP7555g7dy5SUlJQX1+PkydPYtasWWhsbGw1hDAqKgojRozAmTNnsHDhQsTHx0OhUGDu3LmYO3cufH198e677+Lxxx/Hs88+ixUrVmDYsGHQaDTIy8vDlStXUF9fj927dyM42D5M7IEHHsDWrVuxdetWTJ8+HVOnToXNZsOBAweQmpoKmUwmDVkk5zVpUiJi4yrxp60XcVlhH1qVrwvEr/KsePC9DXj4oZlQe3EGNepbhrIyvLfuELbq4gGNY/8YWzl+ev8Y+Pn07FlbIqL+xh4qckkTJkzAoUOH8OKLL2Ls2LG4ePEiNm7ciBMnTiAyMhLPPvssDhw4IA1v8/T0xBdffIGf/exn8Pf3x86dO5GRkYHRo0fjr3/9K1555ZV+y37PPffgv//9L+RyObZv346CggLMmzcPGzduREJCQqvj33jjDTz11FPw9PTE/v37kZeXh5/85CdYuXJluz/jX//6F+bPn4+bN2/i888/x8cff4ysrCzp9QkTJmD37t340Y9+BK1Wi0OHDmH//v2oqanB7Nmz8c9//hPx8Y51i2QyGf72t7/h+eefR0BAAPbs2YPz58/je9/7Hv79739zwgoXEhLkiz98OwX3+9ZC9lWvlFUmx2fqYfjFJyeRk3lGcEJyZ9mHjuGnX16yF1NfUVotWBLSiP/37UkspojIJclsTjTOIzs7G2vXrsXly5dhNpsRGRmJefPmIS0trUvXqa6uxvHjx3H16lVcvXoVt27dgtVqxTPPPIMpU6a0eU5JSQlOnDiB06dPIy8vD+Xl5dDpdIiNjcW8efN6faHVgoKCOw6xsdlsKCsrQ0BAQK9/YBU2y98AtmzZMnz22Wf47LPPMHnyZNFxOtTX90df3tvUeeev5uH/DuWjSOkl7VNazXhUeRv3PjATSo2m1TlBQUEA0K0JUci9dXRvGOvr8Onn+7BWPhhWmeO73GhLFZbNjMWQqKB+y0n9j783qCOi7w+ZTNbmMjhd4TRD/jIyMrBixQrYbDaMGDECXl5eOHfuHFauXImbN2/iscce6/S1Ll26hH/84x9d+vnvvPMOsrOzoVarMXToUMTFxaGoqAhnzpyRhk51JQMRkbMbFReBt6OC8d7G49je6A8AMMuVWGUdhOOrD+GpaTGIHDZEcEpyddfOXMSfj5fghs5xL8lsVtznXYNHFyRDreRMk0Tk2pyioKqtrcXKlSthtVrx3HPPISUlBQBQWVmJl156CZs3b8b48ePbHA7VFl9fX8ydOxdxcXEYMmQI1q9fj/T09A7PCQwMxPTp0zFt2jRotVpp/8mTJ/HGG29g8+bNGDt2LMaMGdP9N0pE5GT0GhV+/NBkpJy+ir+crkS5yr50wEVdOJ7JaMBD5/bhwYWToda2XtSZqCP1tXX4eFMGNplDYNUFS/tDTVV4ZlIYRg5re9kGIiJX4xTPUO3Zswf19fVITk6WiinAXhg1zS62adOmTl8vPj4eTzzxBNLS0hAVFdWpIUXPPPMM5syZ06KYAoCkpCTMnDkTgH2aaiIidzR+bBz+76GRmCYrlvaZ5Up8YgjFso+O49zhTM4ESJ1is9lwbM8R/OSzc9hgCWsxxG+eqgRvfXMcRg6LFpiQiKh3OUVBlZmZCcAxRXRzSUlJUKlUyMrKgtFo7O9oAOwzrQFARUWFkJ9P7uHtt99GXl6e0z8/RQOXt6ceP3t0On45uBH+xmppf54mAL++4YE//2crytpYE42oScGFC3j939vxhwI/lKp9pP1hjRV4eQTw44enQf//27vz+Kiq+//jr5nMZJkkJCxJICAEAhpZZI8ICi4gIqhoZdH+bLGlPtAqIorlUazVh1osFLBqq/2qKAU3rGirUlQQkIAsCiiSELNCMIEsZBJISDLb7w+aSEyCYXInk4H38x/h3jtnztWPd87nnns/J0yznSJybmkTCdWhQ4cA6NWrV4N9FouF7t2743A46q2L05pq192Jior6iSNFRALfZSMH8vzNSUzke0wed932DaG9+Pm76az5v1W4Kiv82ENpa5xldlY+9xr/b+33fGFLqNtucTuZYsnnmWmXMGhIkv86KCLiQ35/h6qyspKKilM/zB06dGj0mA4dOpCVlUVxcTEJCQmt2DuoqKioe//q9AVXf8rcuXMbbAsODubpp58GTr2z9VM8Hg+lpaVYLBbDK6HVtldbzU3kdL6OD4/Hg9lsJiYmRlX+2qiYmBgWzEnkhr1pLF6fQZbl1LpV5dYI/lIZQZ9VO7h3QBTDJl2Lydwm7s2JH3hcTj5/9yP+nukkz5ZYb1TRz1XC7yZeQu+kq/zXQWkTan9Laqu5iZzuXIgPv4+mq6qq6v4c0kiJ3tO3n35sa3nppZcoLy+nT58+JCcnt/r3i4j40yWDLmb5gIt46/3PWX7QTbXZCkBGeDz3Z8PlT6/it9cPosegS/zcU2lt323bwXObMvnK1h1sP2wPd1Vx98XhTJ5wE2bdMBGR84AhCdWSJUvIy8s7q8/ce++99O7d24iv95n333+fbdu2ERERwezZs8/qTvrSpUvPuL+4uLhZ61C53W6cTqfWoZJW1RrrULndboqKijRDFSCuG92PwaUn+OdnaaRU/bBuVYqtJzs+K+aGDSu49fpLCe/U0Y+9lNZQll/Amx/v5mNrAm7bD8UlTB4P4yIruf2afrSPCKWkuNiPvZS2xN/rDEnb5u/4aDPrUBUVFZ31+03V1dUA9arqVVdXY7PZmnWsr23atIk333yTkJAQ5s+fT1xcXKt9t4hIWxTXPoJFs67nq29zeObjb8kKigbAYbayhl5s+CCH6e2/YeyEUQSHqPDAuaaq4iRr127lX5UdqQiuvz5ZX08pc28cSt/Erho0i8h5x5CEqva9IG/YbDZsNhuVlZUcO3as0YTq2LFjQPPeOzLCrl27ePHFFwkKCuKhhx7iwgsvbJXvFREJBEP79+QvseFs3LafVZlVHLOcWruqLDiCf1REsOaNPUyJc3D1NcOxNvEotwSOqooKPl7/JWtKw7Fbu9UbOcQ5ypnRP5rLho4gNja26UZERM5hfn+HCk6VJU9LSyM7O5tu3brV2+d0Ojl06BBWq5X4+Hif92X//v0888wzAMyePVsL+YqINMJsMnHNqP6MHOZgzboveb8sgpqgU+9XFQVH8fdSePf13UzpUMlVYy/FEhHh5x7L2aq22/n0053860R7SoPjwPrDPpuriikxNUwaO4Rga5sYSoiI+E2bKM00ZMgQALZv395g3+7du3E4HPTv35/gYN8+QpKdnc2iRYtwOp3MmjWr0XWxRETkB2EhVn5+02X8fUI3ruV7gtyuun1HQ9rzfEVX7n37azas/ghn6TE/9lSaq6boKGvf+IC716TzUk13SoN/eGcu2OXghqACXrgpkVsmJCuZEhGhjSRU11xzDWFhYXz55Zfs2LGjbntZWRmrVq0CYNKkSQ0+N2fOHObMmVP3SGBL5Ofn86c//YmTJ08yY8YMrrzyyha3KSJyvoiJac9vf34Nfx8dzVgKMHt+SKwKQjvyrCORu9cc4N8r/0NF/vd+7Kk0pSwnh9Wv/Zu7PjjIPzx9KAn5Ye1Fq9vBRHMBL46LY+b0q4iOijxDSyIi5xeT56dKzbWS7du3s2zZMgD69u1LZGQk+/bto6KiggkTJnDnnXc2+MzUqVMBeP755xs8u71gwYK6Px85coTjx4/TuXNnIiNP/Qj07NmTmTNn1h3z8MMPk5ubS7t27Rg8eHCjfezatSuTJ09u0XnWKigoaFaVv5KSEjp27Kgqf9KqWqPKn69iW3yrudWY8guKWL05nc3ODrhN9e/d2ZxVjHMfZtLIJGIvatvVXs8H3+9L5T87s/nM2p2aoPpPgljcTsaF2rn1qgF06nTmxe39XalL2i7FhpyJv+OjzVT5M8KIESN4/PHHWbNmDRkZGTidTrp27cr48eO56qqzXxQwIyOjwbYjR45w5MgRAKxWa719tYsLl5eXs3nz5kbb7Nu3r2EJlXiva9euP3nMlClT6t6FE5HWF98lhjnTY7i1sIzVnx9gS1VkXWJVaQnl3/Tmg13VjEpZx40D4+kzdICS61bkdrlI3fE1/04tYlfoBXhC6ye2VreTqyIqmXJVP2Kjw/3USxGRwNBmZqjON5qh8l5tQjVlypQmj0lOTub2229vrS6dtSVLlrB06VKWLl3KtGnT/N2dBjRDJU3x9k5iYdlJPty8j0/tIVQGNaz817O6mHExbkaP7EdkR61l5SulBUfZtP0A60utHA7p0GB/O9dJJnRyMWHMANqHn12FRn/fZZa2S7EhZ+Lv+DinZqhEzpZmoEQCR2xUGL+6MZnp1Q7Wf/4NH3zvotDarm5/Tkgn/q8cXvsonxGOXYxLCKf/ZUMw2zQ70lLO42Xs2bqX9YdPsivkAlzmOPhRrtTVYeemhDDGjLqEUGuQfzoqIhKglFCJiEirsYVYuXHcUK53udmxYz//TreTHhxTt78myMrnQQl8fgQ6v7GXK02FXNG3C12HDcFk1WLBzeU+WUnurj1s/e4oG83xlITEQVjD4y5xFHJj/1iGDr0Us2aLRUS8okf+/ESP/Hmv9pG/779vXqWwJ598khdeeIEbbriBF198sd6+4uJixo4dy7Fjx1izZg3Dhg0D4OjRo7z77rts2LCB3NxcSkpKiI6OZtiwYdx7770MGjSo0e+qrKzklVde4cMPPyQnJwePx0PXrl0ZPXo0d911F926dePSSy/l8OHDjX7+nXfeYeTIkc38N+E7euRPmuKLRzMOpmexfs9BNlVHU25puLg7QPfKo4wMOc7Ift24YGB/zBbdD/wxt6OGrK/2se1AAV8421MQ1vijkx0cJ7jKdpyxw3oT3+sCw77f34/tSNul2JAz8Xd86JG/85jH46HC4fb68xb3qUGs0+n6iSONFW41t/oA+uGHH2bLli188MEHXHPNNfXevXrwwQcpKipi7ty5dckUwMcff8xTTz1FQkICSUlJREREkJuby3//+1/Wr1/PihUrGDNmTL3vOXr0KNOnT+e7774jOjqaUaNGYbFYyM3NZfny5fTr149p06YxceJEtmzZQmpqKsOHDychIaGujR9XqxQ5H/S4KJFfX5TIHQ4Xu746wKeZdvbSAc9p14pDtjgOEcdbB6Dr3h1cFlZBct/uJPbvXXcD4HzkqHHw3Tfp7DjwPV/UtKMwpD1Ye9dbhBcgyONimLmUcRfFMnjQECxBbWLVFBGRc4JmqPykpTNUJ2pc/PydhpUM27rXp/QhIrhlz+ef7QwVQGZmJuPHj8dqtfLpp59ywQUXsGLFCn7/+98zZMgQ3nvvvXqDsrS0NDweD3379q3XzqZNm7jzzjuJj48nJSWl3n+XadOmkZKSwuTJk1m8eDE22w932rOzs3G73fTufaqSlopSaIYqULXWncTC8pNs2vkdXxRUkW1uulx3uPMkA1zFDIyCQYlxdE7qgzm0kWfbzhHuygoOp2WwN7uQr8vNfGuNoaqRIh+1LvTYGdUtgjHD+5x1kYmz5e+7zNJ2KTbkTPwdH5qhkvPamcqnv/LKK1x33XV1f+/duzd/+MMfWLBgAbNnz+bpp5/miSeeIDw8nOeee67BHe6LL7640XavvPJKJk2axJo1azhw4EDdcXv27CElJYXY2NgGyRRAr169vD1NkfNSbLswpo4dyFQgv+Q4X+xKZ9vRGjIt9SvTVVjC2G65gO1VwH6I/WoPlzgKSYoOok/3WLr164MlqmE1u0DhKCnmYGoGmYeKSDvu4ZvgzhwLiQK6Q2jD400eN0muUkZ2DWPE8IuJjUpq9T6LiJxvlFBJwDpT2fTGkq0ZM2bw2WefsWHDBm666SZOnjzJkiVL6j1yd7rq6mo2bdrEnj17OHbsGDU1NcCp2SuAnJycuoRqy5YtANx8880NkikRaZn4jpH87Lph/Aw4UljKF7vS2HnUQbq1Ey5z/RnvwtAOrA/twHoXkAOhmXn0qtpNn1AHfeIiSezTndgL4rFY297Pn7PGwZGDh8nIPExm4QkyqoPJCY2lJqg90B4iG/+c1e3gYmcJl3YJ47Lki+nYoW/jB4qIiE+0vV8UaZZwq5nXp/Tx+vP+KkoRbjXuuX1vyqYvWbKEESNGcPz4ccaNG8f06dMbPS4tLY0777yTvLy8Jts6ceJE3Z/z8/MB6NGjx1n3SUSar3Nse26eOJKbgcoTlez/Nou9eXa+rrSSZ4lucHxVUAip4d1IBSgGiiuxbD1Alxo7XTwVdLU4ibeZiI8KIa5TFFExHbF2isVkcLl2j8cDlRXUFBdSVljCkZJy8suqyT/p4XunlQJTOEeCo/+XIMaAJeaMv9A9naUMjHAyqEdHLu7Xh9CwAYb2V0REmk8JVYAymUwtehfJYjn1Waf5/HqF7pNPPqGqqgo49V5VZWVlgxklj8fDrFmzyMvL44477uCOO+6gR48ehIeHYzKZWLhwIc8//3yj78DpfSCR1mOLsDF8xACGjzj195IT1Xyddoj9ecfIPAGHzJG4TQ1v4jjNFvJCO5FHp1MbHPwv2QIOVBPuSCfaWUmUp5roIDfRwRBqMWPFjRU3wbixmjxYcROEBycmnJip8ZhxYMKBmRrMnHR6sNd4KHMFYTcFU2YNp9ISxqmppv9NN1lpUEDix4I8LhI8x+kdaaZ/945cknQB0TaVkBcRaSuUUMl5Izs7m8ceewybzcbo0aNZt24djz32GIsWLap3XGZmJpmZmQwcOJCnn366QTuHDh1qsC0+Ph6A3Nxcn/RdRH5ax4gQrh7eh6uHn/r7SYeL7JwCMrLzySg+SaYjhCPB0T/ZToXVRoXVRr2yN94URDXRYAHdn/yIx028o4w+wdX0iQmnT+9u9OwRR7Cq8omItFlKqOS84HQ6ue+++6isrGTx4sVMnjyZa6+9ltdff52rr766XgELu90O0GjFF7vdzueff95g+xVXXMGf//xn3nvvPebNm0dY2JmrjFmtp25Ju1ytW7Ze5HwSZg2i34Xd6Hdht7ptVRUnKcgvJL/QzvelJ8mvcJFfYyafMI6bG6ny4CPt3FV0NZ0kPthNfISFru1txMdF0zk+hpDQ1uuHiIi0nBIqOS8sWbKEvXv3Mn78eG6//XYAnnvuOSZPnsy8efMYMmRI3RpQPXv2xGw2s3XrVrKzs+sq9FVVVTF//vy6hOt0gwcPZuTIkWzbto2HH36YRYsW1UuqcnJycLlcdWXT4+LiAMjKyvLlaYvIj4SGh9GzTw969mn4vmONy4290oG9tAx7sR27vRz78ZOUVTio9lDvcT4HQTj+96ifpfYxwP/9s/bvIWYTUbZg2rezERUdSXSn9kR3aEd0qBVrkB4PFhE5VyihkoA1Z86cJvd17dqVefPmAbBz507+9re/ERMTw+LFi+uOGTx4MHPmzOEvf/kLc+fOZeXKlZhMJjp16sRtt93G66+/zrhx4xg1ahShoaHs3LkTl8vF1KlTWb16dYPvfPbZZ5k6dSpr1qxh48aNJCcn1y3sm5qaypIlS+oSqjFjxhAaGspLL71Eeno6cXFxmEwmZs2aVXeMiLSu4CAzsZEhxEbGQnctsi0iIs2jhEoC1jvvvNPkvr59+zJv3jyOHz/O7NmzcblcLFmyhI4dO9Y7bvbs2WzcuJGNGzfy6quv8qtf/QqAhQsXkpiYyFtvvcXWrVuJjIzkiiuu4He/+x1vv/12o9/ZpUsX1q5dy0svvcRHH33E5s2bsVgsxMfHM3PmTC6//PK6Yzt37szy5ctZtmwZO3fupKKiAoBbbrlFCZWIiIhIADF5GitVJj5XUFDQaJW403k8HkpKSujYsaPh1eP8VTZdAoOv48OXsS2+5e8V7aXtUmxIUxQbcib+jg+TydToe/NnQ2WDREREREREvKSESkRERERExEtKqERERERERLykhEpERERERMRLSqhERERERES8pIRKRERERETES0qoREREREREvKSEqg0zmUyYTCbcbre/uyJiKJfLVRffIiIiIoFMCVUbFxwcTFVVlb+7IWKo6upqgoOD/d0NERERkRaz+LsDcmZhYWGUlZUBEBISQlBQkCHtejyeev8UOZ2v4sPlclFdXU1VVRVRUVGGti0iIiLiD0qo2jiLxUJUVBQnT56krKzMsAGu2XxqclKPE0pjfBUfJpOJ4OBgoqKisFh0+REREZHApxFNALBYLERGRgLGzRjExMQAUFRUZEh7cm7xVXzonSkRERE51yihCjBGDUhr29EAVxqj+BARERFpHhWlEBERERER8ZISKhERERERES8poRIREREREfGSEioREREREREvKaESERERERHxkqr8+UlbqZ7WVvohbZPiQ5qi2JCmKDakKYoNORN/xYcR32vyGLWwkYiIiIiIyHlGj/yJiIiIiIh4SQnVeWr+/PnMnz/f392QNkrxIU1RbEhTFBvSFMWGnMm5EB96h+o8VVNT4+8uSBum+JCmKDakKYoNaYpiQ87kXIgPzVCJiIiIiIh4SQmViIiIiIiIl5RQiYiIiIiIeEkJlYiIiIiIiJe0DpWIiIiIiIiXNEMlIiIiIiLiJSVUIiIiIiIiXlJCJSIiIiIi4iUlVCIiIiIiIl5SQiUiIiIiIuIlJVQiIiIiIiJeUkIlIiIiIiLiJSVUIiIiIiIiXrL4uwNijJqaGt5//322bt1KcXExERERDBw4kGnTptGxY8ezaquiooJ33nmHnTt3YrfbiY6OZvjw4UydOpXw8HAfnYH4ihGxUVFRwZ49e/jqq6/Izc2luLgYk8lEt27duPzyy7n22muxWHQ5CURGXjtOV1BQwEMPPYTD4WDgwIEsWLDAwF5LazA6No4cOcL777/Pvn37sNvthIaG0qVLF5KTk7nxxht9cAbiK0bGxt69e1m7di1ZWVlUVlYSHh5O7969mThxIgMGDPDRGYgvZGdn880335CZmUlGRgalpaVYrVZef/11r9oLpPGoyePxePzdCWmZmpoannjiCdLT02nfvj1JSUkUFRWRmZlJu3btePLJJ+ncuXOz2jp+/DiPPPIIBQUFxMXF0atXLw4fPkxeXh6dO3fmqaeeIjIy0sdnJEYxKjbeeust1qxZg8lkomfPnnTu3Jny8nLS09NxOBwkJSWxYMECQkJCWuGsxChGXjt+7PHHHyc1NRWPx6OEKgAZHRs7d+7kr3/9K06nk4SEBLp06cKJEyc4dOgQISEhPPfccz48GzGSkbHx4Ycf8s9//hOTycRFF11Ehw4dOHr0KFlZWQDMnDmTa6+91penIwZatGgRX375Zb1t3iZUgTYe1S3lc8B7771Heno6F154IY888gihoaHADxeqF154gccff7xZba1YsYKCggKSk5N54IEHCAoKAmD58uWsW7eOFStWcO+99/rsXMRYRsVGaGgoN998M+PHj6dDhw512wsKCnjiiSc4cOAA7777LrfffrvPzkWMZ+S143SfffYZ+/fvZ+zYsaxfv97obksrMDI2cnNzeeaZZwgLC2PevHkkJSXV7XO73eTk5PjkHMQ3jIqN8vJy3njjDSwWC48++mi9uNi+fTvLli1j5cqVjB49uu47pG278MILSUhIIDExkcTERO666y6v2wq08ajeoQpwTqeTdevWAfDrX/+63kVn0qRJ9OjRg7S0NLKzs3+yLbvdzpYtWwgKCmLmzJl1wQtwxx130K5dO1JSUrDb7YafhxjPyNiYPHkyt912W71kCqBLly51SdTWrVsN7L34mpHxcbqysjJWrlzJgAEDGDVqlKF9ltZhdGy8+uqrOJ1O7rnnnnqDZgCz2UxiYqJxnRefMjI2MjIycDqd9O/fv0FcjBgxgu7du1NdXc3hw4eNPQnxmcmTJzN16lSGDh1KdHS01+0E4nhUCVWAO3DgABUVFcTFxdGzZ88G+y+99FKABlOwjdmzZw8ej4e+ffs2+B/BarUydOhQ3G43e/fuNaLr4mNGxsaZJCQkAFBaWtqidqR1+So+Xn31VWpqavjNb35jSD+l9RkZG4cPHyYtLY0uXbowdOhQw/sqrcvI2LBarc36zoiIiLPrpAS8QByPKqEKcAcPHgRo9MIG0KtXr3rHtaSt2u25ubln203xAyNj40yOHj0K0KK7UdL6fBEfu3fvZtu2bdx8881ev3sl/mdkbHz77bcAXHLJJdTU1LBp0yaWL1/O8uXL2bBhA5WVlQb1WlqDkbGRmJiIzWbj22+/5cCBA/X27dixg0OHDnHRRRfpWnIeCsTxqN6hCnDFxcUATVbVqX1Eq/a45rT148e6atV+R3PaEv8zMjbOZO3atQAMGzasRe1I6zI6PqqqqnjllVeIj49n8uTJhvRR/MPI2MjLywMgODiYhx9+mPz8/Hr733jjDR588EH69u3bki5LKzEyNsLDw5k1axbPPvssf/zjH+uKUhQWFpKVlcWgQYO45557jOu8BIxAHI8qoQpwVVVVAE1WV6t9vrn2uJa0Vbu9urr6rPsprc/I2GjKJ598wr59+wgPD9cgOsAYHR9vvfUWRUVFPProoyqhH+CMjI2Kigrg1I2X8PBwHnroIfr374/dbudf//oXKSkpLF68mKVLl9K+fXuDzkB8xejrxogRI4iIiGDZsmX1ZqmioqLo169fm6riJq0nEMej+tULcD9V9f5squLXHmsymVrUJ2kbjIyNxqSmpvLaa69hMpm4++67m7yTJG2TkfGRlZXFunXrGD16NP37929p18TPjIwNt9sNgMvl4r777mPgwIEA2Gw2Zs+eTUFBAVlZWXz88cdMnz7d+05LqzD6d+WDDz5g1apVdWsLxcbGUlhYyNtvv82qVavIyMjgwQcfbEmXJQAF4nhU71AFuLCwMKDpLL12e3NKjta21dSdpdq2tNZQYDAyNn7s4MGDLF68GKfTyYwZM0hOTva+o+IXRsWHy+XiH//4BzabjV/84hfGdlL8wshrR+0xHTp0qEumTnfVVVcBsH//fq/6Kq3LyNhITU1l5cqVJCQkMHfuXLp3705oaCjdu3fnwQcfpGfPnuzYsYOvv/7auBOQgBCI41HNUAW4Tp06AVBSUtLo/mPHjtU7rjlt1X7mx2q/ozltif8ZGRunO3LkCE899RQVFRVMmTKFCRMmtKyj4hdGxUdJSQm5ublER0ezdOnSevtqH/fKzMzkscceIzQ0lPnz57e06+JjRl47YmNjAYiJiWl0f+328vLys+6ntD4jY2Pz5s3AqcqAZnP9+/tms5nk5GRycnLYv39/o8m4nLsCcTyqhCrA9ejRA6DJhRFr14KoPa4lbdVub05b4n9GxkatY8eO8eSTT2K327n++uuZMmVKyzsqfmF0fNjt9ibXBKmoqCA1NRWbzXb2HZVWZ2Rs1C6rcOLEiUb3Hz9+HPBuplxan5GxUTtYrp2N+LHa7U3Fjpy7AnE8qoQqwCUlJWGz2Th69Cg5OTkNSkzu2LEDgCFDhvxkW4MGDcJkMpGWlkZZWRlRUVF1+xwOB1999RUmk4nBgwcbexLiE0bGBpz6UXvqqacoLCzkyiuv5Je//KXhfZbWY1R8xMbGsnr16kb37d+/n8cff5yBAweyYMECYzouPmfktWPAgAGEhIRw5MgRiouLG9xRTk1NBZoujyxti5GxUTvGyMrKanR/7fbaWU45fwTieFTvUAU4i8XCddddB8Dy5cvrPW/64YcfcvDgQZKSkujdu3fd9nXr1jFnzhzeeOONem21b9+eUaNG4XQ6efnll3G5XHX7Vq1aRXl5OZdffrnWGwoQRsZGdXU1CxcuJC8vj8suu4xZs2YF1Mui0pCR8SHnFiNjIyQkhAkTJuByuXj55ZfrtbV37142b96MyWRi7NixPj4rMYKRsVH77m1KSkqDhYB37dpFSkoKJpNJ7+iew86l8ahmqM4Bt9xyC/v27SM9PZ3777+fpKQkiouLycjIIDIyssE6DuXl5eTn51NaWtqgrRkzZpCRkcGOHTuYM2cOiYmJ5OXlkZeXR1xcnGYlAoxRsfHmm2+SkZGB2WwmKCiIF154odHv++1vf+uzcxHjGXntkHOLkbFx6623kpaWxu7du7n//vvp3bs35eXlfPfdd3g8HqZPn15vAC5tm1GxMXz4cEaMGMH27dtZtGgRiYmJxMTEUFRUVDc7NX36dOLj41vt3KRldu/ezbvvvltvm9PprPeEws9+9rO6GcxzaTyqhOocEBwczB//+Efee+89UlJS2LVrF+Hh4YwZM4Zp06ad1Ut77dq1Y+HChaxevZpdu3axc+dOoqKiuO6665g6dSoRERE+PBMxmlGxUVtcwO12k5KS0uRxSqgCi5HXDjm3GBkbtW395z//YcuWLezduxer1Ur//v2ZOHFisx87lrbBqNgwmUw88MADbNy4kc2bN3Po0CFyc3Ox2WwMHjyYCRMmMGjQIN+ejBiqvLycjIyMets8Hk+9bc0tQBNo41GTp6WL0YiIiIiIiJyn9A6ViIiIiIiIl5RQiYiIiIiIeEkJlYiIiIiIiJeUUImIiIiIiHhJCZWIiIiIiIiXlFCJiIiIiIh4SQmViIiIiIiIl5RQiYiIiIiIeEkJlYiIiIiIiJeUUImIiIiIiHhJCZWIiIiIiIiXlFCJiIiIiIh4SQmViIiIiIiIl5RQiYiIiIiIeEkJlYiIiIiIiJeUUImIiIiIiHhJCZWIiIiIiIiX/j/aPPqPF5VNvwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 960x720 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "plt.style.use('ggplot')\n",
    "plt.rcParams['figure.dpi'] = 150\n",
    "\n",
    "plt.plot(x[1:-1], u)\n",
    "\n",
    "xe = np.linspace(0, 1, 101)\n",
    "plt.plot(xe, -0.5*xe + 0.5*xe**2)\n",
    "plt.legend([\"Computed\", \"Exact\"])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}